Session Three

Contouring gridded data sets

GMT comes with several utilities that can create gridded data sets; we will discuss two such modules later this session. The data sets needed for this tutorial are obtained via the Internet as they are needed. Here, we will use grdcut to obtain and extract a GMT-ready grid that we will next use for contouring:

G = grdcut("@earth_relief_05m", region=(-66,-60,30,35), verbose=true)

Here we use the file extension .nc instead of the generic .grd to indicate that this is a netCDF file. It is good form, but not essential, to use .nc for netCDF grids. Using that extension will help other programs installed on your system to recognize these files and might give it an identifiable icon in your file browser. Learn about other programs that read netCDF files at the netCDF website. You can also obtain tut_bathy.nc from the GMT cache server as we are doing below. Feel free to open it in any other program and compare results with GMT.

We first use the GMT module grdinfo to see what’s in this file:

grdinfo(G)

The file contains bathymetry for the Bermuda region and has depth values from -5475 to -89 meters. We want to make a contour map of this data; this is a job for grdcontour. As with previous plot commands we need to set up the map projection with proj. Here, however, we do not have to specify the region since that is by default assumed to be the extent of the grid file. To generate any plot we will in addition need to supply information about which contours to draw. Unfortunately, grdcontour is a complicated module with too many options. We put a positive spin on this situation by touting its flexibility. Here are the most useful options:

Option	Purpose
annot=annot_int	Annotation interval and attributes
cont=cont_int	Contour interval
labelsgap	Controls placement of contour annotations
range=(low,high)	Only draw contours within the low to high range
cut=cut	Do not draw contours with fewer than cut points
smooth=factor	Resample contours smooth times per grid cell increment
ticks=(local_high=true, local_low=true, gap=gap, closed=true, labels=lab)	Draw tick-marks in downhill
	direction for innermost closed contours. Add tick spacing
	and length, and characters to plot at the center of closed contours
pen=(annot=true, contour=true, pen=pen, colored=true, cline=true, ctext=true)	Set contour and annotation pens
scale=(factor=factor, shift=shift, periodic=true)	Subtract offset and multiply data by factor prior to processing

We will first make a plain contour map using 1 km as annotation interval and 250 m as contour interval. We choose a 15-cm-wide Mercator plot and annotate the borders every 2°:

grdcontour(“@earth_relief_05m”, region=(-66,-60,30,35), proj=:Mercator, figsize=15, cont=250, annot=1000, show=true)

Your plot should look like our example 11 below

Result of GMT Tutorial example 11

Exercises:

1. Add smoothing with smooth=4.

2. Try tick all highs and lows with ticks.

3. Skip small features with cut=10.

4. Override region using region=(-70,-60,25,35)

5. Try another region that clips our data domain.

6. Scale data to km and use the km unit in the annotations.

Gridding of arbitrarily spaced data

Except in the situation above when a grid file is available, we must convert our data to the right format readable by GMT before we can make contour plots and color-coded images. We distinguish between two scenarios:

1. The (x, y, z) data are available on a regular lattice grid.

2. The (x, y, z) data are distributed unevenly in the plane.

The former situation may require a simple reformatting (using xyz2grd), while the latter must be interpolated onto a regular lattice; this process is known as gridding. GMT supports three different approaches to gridding; here, we will briefly discuss the two most common techniques.

All GMT gridding modules have in common the requirement that the user must specify the grid domain and output filename:

Option	Purpose
region=(xmin,xmax/ymin,ymax)	The desired grid extent
inc=(xinc [,yinc])	The grid spacing (append m or s for minutes or seconds of arc)
outgrid=gridfile	The output grid filename

Nearest neighbor gridding

Search geometry for nearneighbor.

The GMT module nearneighbor implements a simple “nearest neighbor” averaging operation. It is the preferred way to grid data when the data density is high. nearneighbor is a local procedure which means it will only consider the control data that is close to the desired output grid node. Only data points inside a specified search radius will be used, and we may also impose the condition that each of the n sectors must have at least one data point in order to assign the nodal value. The nodal value is computed as a weighted average of the nearest data point per sector inside the search radius, with each point weighted according to its distance from the node. The most important switches are listed below.

Option	Purpose
search_radius=val	Sets search radius. Append unit for radius in that unit [Default is x-units]
empty=*val*	Assign this value to unconstrained nodes [Default is NaN]
sectors=*n*	Sector search, indicate number of sectors [Default is 4]
weights=true	Read relative weights from the 4th column of input data

We will grid the data in the file tut_ship.xyz which contains ship observations of bathymetry off Baja California. We obtain the file via the cache server as before. We desire to make a 5’ by 5’ grid. Running gmt info on @tut_ship.xyz yields:

tut_ship.xyz: N = 82970     <245/254.705>   <20/29.99131>   <-7708/-9>

so we choose the region accordingly, and get a view of the contour map using

G = nearneighbor("@tut_ship.xyz", region=(245,255,20,30), inc="5m", search_radius="40k");
grdcontour(G, proj=:Mercator, figsize=15, cont=250, annot=1000, show=true)

Your plot should look like our example 12 below

Result of GMT Tutorial example 12

Since the grid ship.nc is stored in netCDF format that is supported by a host of other modules, you can try one of those as well on the same grid.

Exercises:

1. Try using a 100 km search radius and a 10 minute grid spacing.

Gridding with Splines in Tension

As an alternative, we may use a global procedure to grid our data. This approach, implemented in the module surface, represents an improvement over standard minimum curvature algorithms by allowing users to introduce some tension into the surface. Physically, we are trying to force a thin elastic plate to go through all our data points; the values of this surface at the grid points become the gridded data. Mathematically, we want to find the function z(x, y) that satisfies the following equation away from data constraints:

\[(1-t)\nabla ^2 z - t \nabla z = 0,\]

where t is the “tension” in the 0-1 range. Basically, for zero tension we obtain the minimum curvature solution, while as tension goes toward unity we approach a harmonic solution (which is linear in cross-section). The theory behind all this is quite involved and we do not have the time to explain it all here, please see Smith and Wessel [1990] for details. Some of the most important switches for this module are indicated below.

Option	Purpose
aspect_ratio=val	Sets aspect ratio for anisotropic grids.
convergence=limit	Sets convergence limit. Default is 1/1000 of data range.
tension=val	Sets the tension [Default is 0]

Preprocessing

The surface module assumes that the data have been preprocessed to eliminate aliasing, hence we must ensure that this step is completed prior to gridding. GMT comes with three preprocessors, called blockmean, blockmedian, and blockmode. The first averages values inside the grid-spacing boxes, the second returns median values, wile the latter returns modal values. As a rule of thumb, we use means for most smooth data (such as potential fields) and medians (or modes) for rough, non-Gaussian data (such as topography). In addition to the required region and inc switches, these preprocessors all take the same options shown below:

Option	Purpose
reg=true	Choose pixel node registration [Default is gridline]

With respect to our ship data we preprocess it using the median method:

D = blockmedian("@tut_ship.xyz", region=(245,255,20,30), inc="5m", verbose=true);

The output data can now be used with surface:

G = surface(D, region=(245,255,20,30), inc="5m", verbose=true);

If you rerun grdcontour on the new grid file (try it!) you will notice a big difference compared to the grid made by nearneighbor: since surface is a global method it will evaluate the solution at all nodes, even if there are no data constraints. There are numerous options available to us at this point:

1. We can reset all nodes too far from a data constraint to the NaN value.

2. We can pour white paint over those regions where contours are unreliable.

3. We can plot the landmass which will cover most (but not all) of the unconstrained areas.

4. We can set up a clip path so that only the contours in the constrained region will show.

Here we have only time to explore the latter approach. The mask module can read the same preprocessed data and set up a contour mask based on the data distribution. Once the clip path is activated we can contour the final grid; we finally deactivate the clipping with a second call to mask. Note also that since we are appending layers to the figure, second and on commands use the bang (!) form. Here’s the recipe:

D = blockmedian("@tut_ship.xyz", region=(245,255,20,30), inc="5m", verbose=true);
G = surface(D, region=(245,255,20,30), inc="5m", verbose=true);
mask(D, region=(245,255,20,30), inc="5m", figsize=15)
grdcontour!(G, cont=250, annot=1000)
mask!(end_clip_path=true, show=true)

Your plot should look like our example 13 below

Result of GMT Tutorial example 13

Exercises:

1. Add the continents using any color you want.

2. Color the clip path light gray (use fill in the first mask call).