Data files and initializing a project

Data on disk

Before using the spacegrids module for data analysis, some minimal organizing of your Netcdf data files inside directories is helpful. The module looks for projects inside the ~/PROJECTS directory tree. Create ~/PROJECTS and mark a subdirectory of ~/PROJECTS, say test_project as a project directory simply by including a text file named projname containing the name of the project (e.g. echo foo > projname for a project named foo). Passing the nonick = True argument to sg.info and sg.Project removes the need for the projname file (see below). This is recommended for initial use. Subdirectories, and Netcdf files placed alongside them, in such a project directory (test_project) are interpreted as experiment output. So both a directory and a Netcdf file can be associated with an experiment. The case where a subdirectory represents an experiment allows the added advantage of being able to store the configuration files that were used by the numerical (Fortran or C) experiments alongside their Netcdf output files. There is a natural framework in spacegrids for parsing these configuration files and adding the parameters to exper objects, allowing the subsequent creation of new coord and axis objects across experiments from these parameters (e.g. CO2 concentrations, specified in the configuration file, might vary across experiments: maybe you would like to see how the average air temp and other quantities depend on it, casting pCO2 as a Coord). However, basic productive use of spacegrids is much simpler than this.

We will analyze some climate model data using a ~/PROJECTS tree containing two numerical model output files (experiments) and one file containing observational data (also captured in an experiment object), in this case each placed within their own experiment directory. This example project tree can be obtained by: wget http://web.science.unsw.edu.au/~wsijp/code/examples.tar.gz. (Unpacking is done inside the $HOME directory using tar xvzf examples.tar.gz, unless a different root path is set.) This provides a project called test_project (note the text file called projname). The UVic climate model output consists of equilibrium fields where the model has been run with two geographic configurations. One where the Drake Passage has been closed by introducing an artificial land bridge (the DPC subdirectory), and one where the passage remains open (the DPO subdirectory). In addition, there is some observational data, the file Lev.cdf placed alongside the two experiment directories, to compare the model results to.

Getting started

With the directory tree in place, we can start analyzing the data.

Overview of the available data

The info function returns a dictionary of all the available project names and their paths. Here, we have one project (named my_project). Make sure there is at least one project listed when following these examples.

Using ls to examine the project directory content for test_project, we see two experiment directories and one Netcdf file. The sg module will interpret all three as experiment objects.

To obtain a summary of the available projects while obtaining the paths, use D = sg.info(). To start a project instance, named P, it is easy to use:

Ignore the harmless warnings. They relate to the interpretation of the axes and dimension names found in the Netcdf files. Note that the Lev.cdf experiment object is named after the file, whereas the other two are named after the directories. Passing the nonick = True option would have been: D = sg.info(nonick=True) and P = sg.Project(D['my_project'],nonick = True). In this case, no projname file would be required, and the project would receive the same name as the directory, unless the name is specified. Upon project initialization, meta data associated with the grids, dimensions and axis directions is collected from the data files and interpreted. The project is now initialized. Subdirectories of the test_project directory correspond to (numerical) experiments and also Exper objects. If Netcdf files had been present in the test_project directory, they would have been interpreted as experiments. In this case, an Exper object (experiments) corresponds not to subdirectory, but a single Netcdf file. This latter arrangement corresponds to the structure of a project saved with P.write(), where all experiments are written to experiment files. If a path to a file (not project directory) is provided to sg.Project, a project is initialised assuming all project is contained in that single file.

The Field class

Fields and their grids

To look at an example of fields (e.g. temperature or humidity), it is best to load one into the project.

The ocean temperature Field is now loaded into memory and associated with the project P. How do you know what variable name to pick? Inspection on the Bash command line of the Netcdf file with ncdump -h usually provides this information. Alternatively, experiments have an ls function: P['DPO'].ls() would yield a long list of available variables to be loaded as fields.

Notice that O_temp has not been loaded for Lev. This is because the Netcdf file inside the Lev directory follows different naming conventions. Inspection of the Netcdf file in the Lev directory reveals that temperature is referred to as otemp. To load it, we use:

Now all three temperature fields are loaded. We access the temperature fields of the DPO experiment and the observations Lev and test resulting field instances:

The shape of a Field corresponds to the shape of the Numpy array it contains and is consistent with the shape of its grid:

The names of the grid components depth, latitude and longitude correspond to the names used in the Netcdf file. For instance, there are different naming conventions for the Lev data:

In fact, the Coord (dim) names in the Lev.cdf file were X, Y, Z, but since these names coincide with the axes names, they have been automatically renamed with the _crd suffix to distinguish them from the axes:

Fields can be saved to Netcdf again, along with all their metadata and coordinates via their write method, e.g. TEMP2.write(). To specify the path for instance: TEMP2.write(path = '/home/me/'). Experiment and coord objects also have save methods. In the case of an experiment, all loaded fields are saved to one file. An entire project can be saved with P.save(), yielding a project file structure (directory containing a projname text file and Netcdf files corresponding to experiments). Generally, it is good to provide a name argument different from the project name, and safeguards exist to avoid overwriting the original project on disk, but no safeguards exist yet to check for multiple instances of the same project name (defined in the projname file): for now, this is up to the user to check. Finally, an arbitrary list of fields can be saved with sg.nugget.

An overview of the experiments and loaded fields corresponding to the project are shown by:

And for an experiment:

*** Warning: whenever the sg module is reloaded (e.g. with reload(sg)), all the above steps and the below steps need to be done again in sequence. Stale objects will lead to strange errors! ***

Slicing, plotting and interpolation

Objects representing general coordinate directions (X,Y axes etc) have been constructed during initialization of project P, and can be accessed as follows:

Bring these axes into the namespace under their own names:

So that we can reference them:

Field objects allow slicing by reference to axis name and using ndarray index values. Here, we use the axis object, followed by a slice object (e.g. TEMP[X,5] or TEMP[Z,:] or TEMP[Y,55:] or TEMP[X,slice(55,None,None))]. We will refer to the use of Ax and Coord objects in __getitem__ as "hybrid slice notation". To plot the surface temperature using sg version of contourf:

The sg module automatically selects sensible contour levels and x,y ticks. Versions older than 1.1.4 lack this functionality. The number of levels can be set using argument num_cont. Alternatively, the Matplotlib levels arguments can be passed. Labels and orientation are extracted from the grid metadata. The field grid is used to extract plotting information, unless another grid is specified in the grid argument.

Slicing can also be done via indices using the normal ndarray notation:

Another note of warning on stale objects. Slicing using Ax objects is a common case where stale objects after a module reload lead to strange errors. If we were to reload the module in the above case, and bring the X,Y,Z,T Ax objects into the namespace again, but use the now stale old TEMP field, the above slicing with Z will lead to a perplexing error. All objects will need to be refreshed in this case, which means redoing all the previous steps. In most common cases a module reload is not required, and this type of confusion arises mostly in module code development situations.

The load method (of the Project and Exper classes) also allows memory loading of subsets of the data only by passing the slices argument:

Means are calculated easily by division.

Basin-wide means are calculated using the maskout method to mask out all values outside the domain enclosing the node (a x,y location) argument. maskout uses a floodfill algorithm to mask the basin up to its boundaries. Note the hybrid notation for the node coordinates.

Plotting these two means for the upper 1600m or so:

A calculated field such as mTEMP can be inserted into the project using the insert method. This is generally not needed, but can be useful when the project or experiment is later written to disk. The insert method takes (name,value) pairs. If the value is a field, it is inserted into the vars attribute under name, otherwise into params.

Now, mtemp is part of the DPO experiment:

To see the effect of opening Drake Passage in the climate model:

Data is generally defined on different grids, so interpolation is needed. In our example, the observations are defined on a finer grid than the model data:

To overlay a plot of the model geography on atmospheric fields, we can load the ocean bathymetry field G_kmt, interpolate it on a finer grid with sg.finer_field (to avoid angled contours and show true resolution) and use contour plot over the 0.5 contour to capture the transition from land (value 0) to ocean(values 1 to 19). Here, we overlay this geographic outline over sea level air temperature:

To interpolate a field F in Spacegrids, simply call the regrid method of F on the desired grid gr as F(gr). For instance, to compute the difference between the zonal mean of the model and the observations:

As a shorthand, call the regrid method simply by calling the field itself on the target grid:

The model is pretty good! It is generally a bit too warm (perhaps due to the use of post-industrialized CO2). The resulting grid is defined on that of Lev:

We have seen that it is easy to regrid data. This is also the case for data defined on grids of different dimensions. The ocean exhibits siginificant variations from its mean with longitude, leading to differences in local climate. How would we compare the sea surface temperature to its zonal mean? The following is bound to give an error, as the zonal mean is defined on a different (smaller dimensional) grid than temperature:

The error message provides a hint: regrid one term. regridding mTEMP[Z,0] onto sTEMP.gr expands the right term (by creating equal copies in the X direction, this is called broadcasting) to yield the right shape:

The North Atlantic is warmer than at similar Pacific latitudes, leading to a milder western European climate. A pronounced zonal temperature difference is also apparent in the tropical Pacific (related to the "warm pool").

To compute the meridional streamfunction, we load the velocity in the y direction and apply the calculations in simple notation:

Now that we have the required fields, we can plot them.

There is strong southern sinking in the DP closed case, and a switch to northern sinking in the DP open case. That's why there was warming in the NH.

Again, the Atlantic MOC can be computed by using the maskout of Fields, using V_masked = V[Y,33:].maskout(node = (X,85,Y,30) ) instead of V (see above description for zonal temperature mean). See our special recipy for this.

Field concatenation is done with the concatenate function, similar to Numpy. If we split surface air temperature SAT latitudinally into three pieces:

And re-assemble along the default axis ax = None, sg will guess to use the incomplete (differing) axes:

(This yields the same result as sg.concatenate([SAT1,SAT2,SAT3],ax=Y).) Drawing the zonal mean shows that SAT has been re-assembled correctly:

If multiple files are found in an experiment directory, and they contain the same field, sg will assume they are pieces of the same field and concatenate them upon loading. An example would be files containing different time slices of the same field (e.g. temperature).

Tip: field objects have a squeeze and unsqueeze method: use help to learn more about them. Fields are squeezed upon loading and unsqueezed when saved.

Grids and coord objects

Grid (Gr) objects model the grids on which the field data is defined. They consist of tuples containing Coord objects. Coord objects consist of points along an axis and are read and named from the Netcdf files. To bring these objects into the namespace under their own name:

This allows us to reference the Coord objects built from the Netcdf file by the names extracted from the metadata:

The result of the __getitem__ method here is the ndarray containing the data of the coordinate points. It is important to realize that these coord names (depth, latitude etc) have been read from the Netcdf file and have not been hard coded into Spacegrids.

A shorthand for grid construction from Coord objects is via multiplication:

The latter 3 dimensional grid is associated with ndarray data where the first index represents latitude, the second longitude and the third depth. Field objects contain a grid object that is always congruent with the ndarray holding the field data. The grid objects contained in fields are used to always ensure data alignment, for instance under multiplication.

Coord objects and grid objects are not the same. To construct a one dimensional Gr (grid) object:

Axis alignment

We had:

With data in the ndarray:

where we recognise the length of the depth coord as 19:

Rearrangement of the coordinates:

Data axis alignment is automatic. For instance, multiplying temperature with velocity:

Multiplication alignment works regardless of the order of the coordinates:

Here, V is actually defined on a different grid: the velocity grid:

Why then is the result of TEMP*V defined on the temperature grid? This is because multiplication automatically triggers interpolation of the right multiplicant to the grid of the left multiplicant. Vice versa:

Axis objects

The objects P['DPO'], P['DPC'] and P['Lev'] are experiment objects, and correspond to the directories containing the Netcdf files. Different data ("experiment") files may contain data defined on different grid and use different coordinate naming conventions. A list of the coord objects associated with an experiment object is accessed from the .cstack attribute:

The depth coordinates differ for the model and the observational data.

However, both coord objects point in the same direction: Z. The name ("Z") of this direction has also been read from the Netcdf files. This is encapsulated in the ax class (axis).

Multiplication of a Gr (grid) with an Ax object (representing an axis) picks the Coord component along that axis.

We also saw the Ax object being used to calculate zonal means. In the following example, we load the surface heat flux, demonstrate some operations using Ax objects and compute the poleward heat transport.

Again, ignore the harmless warning: the Lev data contains no heat flux data. The heat flux data is 2 dimensional, as can be seen from the grid:

And looks like this:

Shift fields by reference to the Ax objects:

Note that the grid is updated along with the data (the zero meridian still passes through Greenwich), resulting in a different grid (longitude_rolled appears in place of longitude). This behaviour can be disabled.

Compute the zonal integral (using multiplication notation HF*X) followed by the primitive in the Y direction (using the pipe | notation):

Interesting, the Drake Passage closed configuration has more southward poleward heat transport. We could have used the coord names to compute the PHT:

The axis names X,Y,... are a shorthand for these coord names, and are easier to use.

More grid operations

Operations such as derivatives and integration are methods of the coord class:

Do we get the same result from both shifts? This is a shorthand to access the X coord of HF_shifted2:

The two shifted coord objects are different objects

But contain the same values:

Other derived Coord objects are obtained from slicing:

Similar to the roll function, slicing gives a new grid:

The sum method of a Coord object only sums the ndarray content of a Field:

Similar to Pandas, the NaN (np.nan) values are not counted in these operations. Similarly, the vsum method is the sum weighted with the volumes (lengths, areas) of the grid cells:

Similar for cumsum and vcumsum.

Grid objects also have volume weighted sum methods. These methods are preferred over their Coord counterparts as coordinate cell dimensions may depend on other coordinates (e.g. the distance between longitudinal arcs diminishes towards the poles), and grid objects take care of this.

The volume of the domain of non NaN values inside TEMP is found from:

This can also be achieved with the multiplication notation:

The ones method of Field yields a new Field, where all non- NaN values (of TEMP in this case) are set to 1 and the NaN values remain NaN. Here, the one values constitute the domain of ocean water, as NaN values indicate land (rock).

This allows a calculation of the mean ocean temperature in the model:

This can also be achieved with the division notation:

The Ax class has a derivative method:

Another notation, involving Ax objects, uses the ^ (carat) symbol:

The volume weighted cumulative sum method vcumsum of the Coord class, the primitive integral along that axis, is easily accessed with the pipe | notation and the ax class (here the instance Z). In the latter case the name of the particular coord need not be known.

Operators, divergences and vector fields.

Spacegrids allows representations of vector fields (e.g. velocity fields). We're loading the zonal velocity to examine a 2 dimensional vector fields constructed from U and V.

Earlier, multiplication of two fields yielded a new field. In contrast, multiplying the two fields U and V does not result in a new field here. Instead, we get a container class holding the two fields. This container is the vfield class (vector field).

Why do these fields behave differently under multiplication? Fields have a direction attribute:

The velocity fields U and V have directions X and Y, obviously indicating the direction they are pointing in: they are "directional" fields. In contrast, temperature TEMP is marked as a scaler, this is a "scalar" field. In general, fields of like direction multiply to yield a new field containing the multiplied nparray data. For instance:

This gives a vector field of the velocity components squared:

Vector fields have a direction method that yields a container of the directions of its component fields.

The sg module has its own wrapper to the quiver plot method of Matplotlib that takes vector fields.

The (field) Operator class and its derived classes (e.g. sg.Der, sg.Integ, sg.FieldMul etc.) allows the construction of operators that act on fields and vector fields. E.g. div and curl operators, etc. Ideas for useful operators are shown below. Operators, say K1 and K2, can be multiplied before they are applied. E.g. if K = K1*K2, no evaluation has occurred yet, and evaluation yields: K*F = K1(K2(F)). So K2 acts on F first, then K1 acts on the resulting field. the curl of a 2D vectorfield U*V would then be curl*(U*V). If U*V is defined on 3D space, the curl operator acts on each z-level. For instance, the X-derivative object can be instantiated as follows:

This is equivalent to calling the ddX field operator object on TEMP:

Applying field operators in succession:

Where:

This is the same as successive calls of ddX on the Field:

We could have obtained the Z derivative of the vertical temperature profile above with an operator class:

Another useful field operator is Pick. It picks a component of a vector field:

Here are definitions of some useful field operators:

An advanced tutorial is here.