Interpolate a 3-D cube, 2-D grids or 1-D series from a 3-D data cube or stack of 2-D grids


gmt grdinterpolate cube | grd1 grd2 … -Goutfile [ -D[+xxname][+yyname][+zzname][+vvname][+sscale][+ooffset][+ninvalid][+ttitle][+rremark] ] [ -Etable|line ] [ -Fl|a|c|n[+d1|2] ] [ -Rregion ] [ -Sx/y|pointfile[+hheader] ] [ -T[min/max/]inc[+i|n] |-Tfile|list ] [ -V[level] ] [ -Z[levels] ] [ -bbinary ] [ -dnodata[+ccol] ] [ -eregexp ] [ -fflags ] [ -ggaps ] [ -hheaders ] [ -iflags ] [ -nflags ] [ -oflags ] [ -qflags ] [ -sflags ] [ -:[i|o] ] [ --PAR=value ]

Note: No space is allowed between the option flag and the associated arguments.


grdinterpolate reads a single 3-D netCDF data cube (or a set of 2-D layers) and interpolates along the 3rd dimension for one or more output levels. The data cube must be organized with one or more layers representing the (common) x and y dimensions while the 3rd dimension may represent distance or time; we refer to this dimension as the level. The output layers may be written as a single 3-D cube or as a set of 2-D layers. Alternatively, we interpolate the cube along the level-axis at one or more arbitrary (x/y) coordinates (-S), resulting in a data table with one or more level-series, or we slice the 3-D cube along an arbitrary vertical slice and write that 2-D slice to a grid file (-E).

Required Arguments


Name of a 3-D netCDF data cube to be interpolated. Alternatively, with -Z, you can specify a set of 2-D grid layers instead.


This is the output 3-D data cube file. If -T only selects a single layer then the data cube collapses to a regular 2-D grid file. If outfile contains a C-language format statement for a floating point number (e.g., layer_%6.6f.grd) then we write a series of 2-D grid files which will contain the values for each level [Default is a 3-D data cube, unless only a single layer is implied by -T]. Also see -S for a similar use to write individual level-series tables.

Optional Arguments


Give one or more combinations for values xname, yname, zname (3rd dimension in cube), and dname (data value name) and give the names of those variables and in square bracket their units, e.g., “distance [km]”), cpt to set a default CPT for this grid [turbo] or give - to remove any default CPT already set, scale (to multiply data values after read [normally 1]), offset (to add to data after scaling [normally 0]), invalid (a value to represent missing data [NaN]), title (anything you like), and remark (anything you like). Items not listed will remain untouched. Give a blank name to completely reset a particular string. Use quotes to group texts with more than one word. If any of your text contains plus symbols you need to escape them (place a backslash before each plus-sign) so they are not confused with the option modifiers. Alternatively, you can place the entire double-quoted string inside single quotes. If you have shell variables that contain plus symbols you cannot use single quotes but you can escape the plus symbols in a variable using constructs like ${variable/+/\+}. Note that for geographic grids and cubes (-fg) xname and yname are set automatically. Normally, the data netCDF variable is called “z” (grid) or “cube” (data cube). You can name this netCDF variable via +vvarname.


Specify a crossectinonal profile via a file or from specified line coordinates and modifiers. If a file, it must be contain a single segment with either lon lat or lon lat dist records. These must be equidistant. Alternatively, the format of each line is start/stop, where start or stop are lon/lat (x/y for Cartesian data). You may append +iinc to set the sampling interval; if not given then we default to half the minimum grid interval. If your line starts and ends at the same latitude you can force sampling along the parallel with +p [great circle]. For a line along parallels or meridians you can add +g to report degrees of longitude or latitude instead of great circle distances starting at zero. Append +x to compute distances along a loxodrome (rhumbline) instead of great circle. Instead of two coordinates you can specify an origin and one of +a, +o, or +r. The +a sets the azimuth of a profile of given length starting at the given origin, while +o centers the profile on the origin; both require +l. For circular sampling specify +r to define a circle of given radius centered on the origin; this option requires either +n or +i. The +nnp modifier sets the desired number of points, while +llength gives the total length of the profile. Also note that only one distance unit can be chosen. Giving different units will result in an error. If no units are specified we default to great circle distances in km (if geographic). If working with geographic data you can use -j to control distance calculation mode [Great Circle]. Use -G to set the output grid file name.


Choose from l (Linear), a (Akima spline), c (natural cubic spline), and n (no interpolation: nearest point) [Default is -Fa]. You may change the default interpolant; see GMT_INTERPOLANT in your gmt.conf file. You may optionally evaluate the first or second derivative of the spline by appending +d and specify 1 or 2, respectively.


Specify the region of interest. Using the -R option will select a subsection of the grid. If this subsection exceeds the boundaries of the grid, only the common region will be output. (See full description) (See cookbook information).


Rather that compute gridded output, create time or spatial series through the stacked grids at the given point (x/y) or the list of points in pointfile. If you need a series of points defined by an origin and an end point or similar, you can make such a file first with project. By default we simply sample the cube at each level. Use -T to interpolate the series. The grid level (e.g., depth or time) will be inserted as the third numerical value in the series records, with optional input columns appended, ending with the sampled cube value. Use the optional +h modifier to append header to the trailing text of these input points. On output the trailing text will become the segment header for the series that originate from each point. By default, the table output is written to standard output. Use -G to specify a file name. Alternatively, if you wish each series to be written to its own data file, let the filename in -G have a C-format integer specifier (e.g., %d) and we will use the running point number to create unique file names.

-T[min/max/]inc[+i|n] |-Tfile|list

Make evenly spaced time-steps from min to max by inc [Default uses input times]. For details on array creation, see Generate 1-D Array. Note: If -Z is set and no output times are set with -T we simply rewrite the grid-produced cube as a 3-D data cube file and exit. Also, for -E and -S you may also just give a range via -Tmin/max to limit the layers considered, with no interpolation between the selected layers. If -T is not given and neither -E nor -S are set, then we simply extract all layers within the bounds set by -R.


Select verbosity level [w]. (See full description) (See cookbook information).


Read all 2-D input grids given on the command line and assume they represent the layers in a 3-D cube [Default reads a single 3-D data cube]. Optionally, append levels and assign these to the cube constructed from the grids. The levels may be specified the same way as in -T. If not given then we default to an integer levels array starting at 0.


Toggles between (longitude,latitude) and (latitude,longitude) input/output in -Stable. [Default is (longitude,latitude)].

-birecord[+b|l] (more …)

Select native binary format for primary table input. [Default is 2 input columns].

-borecord[+b|l] (more …)

Select native binary format for table output. [Default is one more than input].

-d[i|o][+ccol]nodata (more …)

Replace input columns that equal nodata with NaN and do the reverse on output.

-e[~]“pattern” | -e[~]/regexp/[i] (more …)

Only accept data records that match the given pattern.

-f[i|o]colinfo (more …)

Specify data types of input and/or output columns.

-gx|y|z|d|X|Y|Dgap[u][+a][+ccol][+n|p] (more …)

Determine data gaps and line breaks.

-h[i|o][n][+c][+d][+msegheader][+rremark][+ttitle] (more …)

Skip or produce header record(s).

-icols[+l][+ddivisor][+sscale|d|k][+ooffset][,][,t[word]] (more …)

Select input columns and transformations (0 is first column, t is trailing text, append word to read one word only).

-n[b|c|l|n][+a][+bBC][+c][+tthreshold] (more …)

Select interpolation mode for grids.

-ocols[,…][,t[word]] (more …)

Select output columns (0 is first column; t is trailing text, append word to write one word only).

-q[i|o][~]rows|limits[+ccol][+a|t|s] (more …)

Select input or output rows or data limit(s) [all].

-s[cols][+a][+r] (more …)

Set handling of NaN records for output.

-^ or just -

Print a short message about the syntax of the command, then exit (Note: on Windows just use -).

-+ or just +

Print an extensive usage (help) message, including the explanation of any module-specific option (but not the GMT common options), then exit.

-? or no arguments

Print a complete usage (help) message, including the explanation of all options, then exit.


Temporarily override a GMT default setting; repeatable. See gmt.conf for parameters.

Generate 1-D Array

We will demonstrate the use of options for creating 1-D arrays via gmtmath. Make an evenly spaced coordinate array from min to max in steps of inc, e.g.:

gmt math -o0 -T3.1/4.2/0.1 T =

Append +b if we should take log2 of min and max, get their nearest integers, build an equidistant log2-array using inc integer increments in log2, then undo the log2 conversion. E.g., -T3/20/1+b will produce this sequence:

gmt math -o0 -T3/20/1+b T =

Append +l if we should take log10 of min and max and build an array where inc can be 1 (every magnitude), 2, (1, 2, 5 times magnitude) or 3 (1-9 times magnitude). E.g., -T7/135/2+l will produce this sequence:

gmt math -o0 -T7/135/2+l T =

For output values less frequently than every magnitude, use a negative integer inc:

gmt math -o0 -T1e-4/1e4/-2+l T =

Append +i if inc is a fractional number and it is cleaner to give its reciprocal value instead. To set up times for a 24-frames per second animation lasting 1 minute, run:

gmt math -o0 -T0/60/24+i T =

Append +n if inc is meant to indicate the number of equidistant coordinates instead. To have exactly 5 equidistant values from 3.44 and 7.82, run:

gmt math -o0 -T3.44/7.82/5+n T =

Alternatively, let inc be a file with output coordinates in the first column, or provide a comma-separated list of specific coordinates, such as the first 6 Fibonacci numbers:

gmt math -o0 -T0,1,1,2,3,5 T =

Note: Should you need to ensure that the coordinates are unique and sorted (in case the file or list are not sorted or have duplicates) then supply the +u modifier.

If you only want a single value then you must append a comma to distinguish the list from the setting of an increment.

If the module allows you to set up an absolute time series, append a valid time unit from the list year, month, day, hour, minute, and second to the given increment; add +t to ensure time column (or use -f). Note: The internal time unit is still controlled independently by TIME_UNIT. The first 7 days of March 2020:

gmt math -o0 -T2020-03-01T/2020-03-07T/1d T =

A few modules allow for +a which will paste the coordinate array to the output table.

Likewise, if the module allows you to set up a spatial distance series (with distances computed from the first two data columns), specify a new increment as inc with a geospatial distance unit from the list degree (arc), minute (arc), second (arc), meter, foot, kilometer, Miles (statute), nautical miles, or survey foot; see -j for calculation mode. To interpolate Cartesian distances instead, you must use the special unit c.

Finally, if you are only providing an increment and will obtain min and max from the data, then it is possible (max - min)/inc is not an integer, as required. If so, then inc will be adjusted to fit the range. Alternatively, append +e to keep inc exact and adjust max instead (keeping min fixed).

File Order

If you provide a series of 2-D files and thus separately assigning the level via -Z, then you must make sure that the order the grids are given on the command line matches the levels you provide via -Z. Unless your files are named in lexical order you must be careful with using wildcards to list all the grids (e.g., *.nc).

Time Coordinates

Time coordinates in netCDF grids, be it the x, y, or z coordinate, will be recognized as such. The variable’s unit attribute is parsed to determine the unit and epoch of the time coordinate in the grid. Values are then converted to the internal time system specified by TIME_UNIT and TIME_EPOCH in the gmt.conf file or on the command line. The default output is relative time in that time system, or absolute time when using the option -f0T, -f1T, or -f2T for x, y, or z coordinate, respectively.

Series creation

The (optional) table-reading and table-producing -S option may require some of the standard common options associated with table i/o, such as -b, -i, o, etc., thus these options are available to grdinterpolate as well. Because the coordinates given via -S are not required to equal the coordinates of the grid nodes, we are resampling each 2-D layer at the given points via grdtrack, hence the availability of the -n option.


To extract a single, new 2-D layer from the 3-D cube for level 3400 using a cubic spline, try:

gmt grdinterpolate -T3400 -Fc

To extract a single, new 2-D layer from the 3-D cube implied by the individual grids layers_*.nc, with individual layer values given via z.txt, for level 3400 using a linear spline, try:

gmt grdinterpolate layers_*.nc -Zz.txt -T3400 -Fl

To resample the 3-D cube for all levels from 1500 to 2500 in steps of 50, using an Akima spline, try:

gmt grdinterpolate -T1500/2500/50 -Fa

The same, but this time write individual 2-D grids per layer, try:

gmt grdinterpolate -T1500/2500/50 -Fa

To extract a time-series through the grids deformation_*.nc at the location (115W, 33N), with the times of each grid provided by the file dates.txt, and append the string “Some like it hot” to the segment header for the series, try:

gmt grdinterpolate deformation_*.nc -Zdates.txt -S115W/33N+h"Some like it hot" > record.txt

To extract a vertical slice of the 3-D grid with seismic velocities that goes through the Hawaii hotspot, selecting cube vs (Isotropic Shear Velocity) and letting the distances be longitude degrees along the parallel, try:

gmt grdinterpolate -E180/20/220/20+i1d+g+p -T25/500/25

See Also

gmt.conf, gmt, grdedit, grdcut, project