Enable writing of 3-D netCDF data cubes from greenspline and grdinterpolate

doc/examples/anim05/anim05.sh

@@ -21,7 +21,7 @@ EOF
 cat << EOF > main.sh
 gmt begin
 	let k=\${MOVIE_FRAME}+1
-	gmt greenspline @Table_5_11.txt -R0/6.5/0/6.5 -I0.05 -Sc -Gt.nc -D1 -Cn\${k} -Emisfit.txt
+	gmt greenspline @Table_5_11.txt -R0/6.5/0/6.5 -I0.05 -Sc -Gt.nc -Z1 -Cn\${k} -Emisfit.txt
 	gmt grdcontour t.nc -C25 -A50 -Baf -BWsNE -JX4i -Gl3.6/6.5/4.05/0.75 -X0.25i -Y0.4i
 	gmt plot misfit.txt -Ct.cpt -Sc0.15c -Wfaint -i0,1,4
 	printf "%2.2d" \$k | gmt text -F+cTR+jTR+f18p -Dj0.1i -Gwhite -W0.25p

doc/examples/ex29/ex29.bat

@@ -16,8 +16,8 @@ gmt begin ex29
 	gmt grdmath -Rg -I4 -r X COSD %a% DIV DUP MUL X SIND %b% DIV DUP MUL ADD Y COSD DUP MUL MUL Y SIND %c% DIV DUP MUL ADD SQRT INV = PROJ_ELLIPSOID.nc
 
 	REM  Do both Parker and Wessel/Becker solutions (tension = 0.9975)
-	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -D4 -Sp -Gmars.nc
-	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -D4 -Sq0.9975 -Gmars2.nc
+	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -Z4 -Sp -Gmars.nc
+	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -Z4 -Sq0.9975 -Gmars2.nc
 	REM Scale to km and remove PROJ_ELLIPSOID
 	gmt grdmath mars.nc 1000 DIV PROJ_ELLIPSOID.nc SUB = mars.nc
 	gmt grdmath mars2.nc 1000 DIV PROJ_ELLIPSOID.nc SUB = mars2.nc

doc/examples/ex29/ex29.sh

@@ -18,8 +18,8 @@ gmt begin ex29
 			SIND $c DIV DUP MUL ADD SQRT INV = PROJ_ELLIPSOID.nc
 
 	#  Do both Parker and Wessel/Becker solutions (tension = 0.9975)
-	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -D4 -Sp -Gmars.nc
-	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -D4 -Sq0.9975 -Gmars2.nc
+	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -Z4 -Sp -Gmars.nc
+	gmt greenspline -RPROJ_ELLIPSOID.nc @mars370.txt -Z4 -Sq0.9975 -Gmars2.nc
 	# Scale to km and remove PROJ_ELLIPSOID
 	gmt grdmath mars.nc  1000 DIV PROJ_ELLIPSOID.nc SUB = mars.nc
 	gmt grdmath mars2.nc 1000 DIV PROJ_ELLIPSOID.nc SUB = mars2.nc

doc/rst/source/api.rst

@@ -1255,7 +1255,7 @@ and pass the ``par`` array with contents as indicated below:
   **GMT_IS_CUBE**.
     Same procedure as for **GMT_IS_GRID** but both ``wesn``, ``inc`` and ``par`` have one extra
     dimension for the depth or time axis.  For non-equidistant layers you need to use
-    ``par[2]`` to sets the number of layers, otherwise ``wesn`` and ``inc`` can set it all.
+    ``par[2]`` to sets the number of layers and use ``inc[2] = 0``, otherwise ``wesn`` and ``inc`` can set it all.
 
   **GMT_IS_DATASET**.
     We allocate an empty :ref:`GMT_DATASET <struct-dataset>` structure consisting of ``par[0]`` tables,
@@ -3795,7 +3795,7 @@ allocated will be *size * n_bands*, where the latter is one of the parameters in
   struct GMT_CUBE {
        struct GMT_GRID_HEADER *header;      /* The full GMT header for the grid */
        float                  *data;        /* Pointer to the float 3-D array */
-       unsigned int           mode;         /* Indicates input was list of 2-D grids rather than a cube */
+       unsigned int           mode;         /* Indicates data originated as a list of 2-D grids rather than a cube */
        double                 z_range[2];   /* Minimum/max z values (complements header->wesn) */
        double                 z_inc;        /* z increment (complements header->inc) (0 if variable z spacing) */
        double                 *x, *y, *z;   /* Arrays of x,y,z coordinates */

doc/rst/source/explain_-D_cap.rst_

@@ -1,7 +1,8 @@
-**-D**\ [**+x**\ *xname*][**+y**\ *yname*][**+z**\ *zname*][**+s**\ *scale*][**+o**\ *offset*][**+n**\ *invalid*][**+t**\ *title*][**+r**\ *remark*]
-    Give one or more combinations for values *xname*, *yname*, *zname* (give the names of those variables
-    and in square bracket their units, e.g., "distance [km]"), *scale* (to multiply grid values after
-    read [normally 1]), *offset* (to add to grid after scaling [normally 0]),
+**-D**\ [**+x**\ *xname*][**+y**\ *yname*][**+z**\ *zname*][**+d**\ *vname*][**+s**\ *scale*][**+o**\ *offset*][**+n**\ *invalid*][**+t**\ *title*][**+r**\ *remark*][**+v**\ *varname*]
+    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]"), *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.
@@ -10,4 +11,6 @@
     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 (**-fg**) *xname* and *yname* are set automatically.
+    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 **+v**\ *varname*.

doc/rst/source/gmt.conf.rst

@@ -949,7 +949,7 @@ PostScript Parameters
         the character set is written as part of the PostScript header.  If
         you need to switch to another character set for a later overlay then
         you must use **--PS_CHAR_ENCODING**\ =\ *encoding* on the command line and
-        not via gmt :doc:`/set`.  Finally, note 6, 8, and 11 do not work with standard fonts.
+        not via gmt :doc:`/gmtset`.  Finally, note 6, 8, and 11 do not work with standard fonts.
 
     **PS_COLOR_MODEL**
         Determines whether PostScript output should use RGB, HSV, CMYK, or

doc/rst/source/grdinfo.rst

@@ -17,9 +17,11 @@ Synopsis
 [ |-E|\ [**x**\|\ **y**][**+l**\|\ **L**\|\ **u**\|\ **U**] ]
 [ |-F| ]
 [ |-I|\ [*dx*\ [/*dy*]\|\ **b**\|\ **i**\|\ **r**] ]
-[ |-L|\ [**0**\|\ **1**\|\ **2**\|\ **p**\|\ **a**] ] [ |-M| ]
+[ |-L|\ [**0**\|\ **1**\|\ **2**\|\ **p**\|\ **a**] ]
+[ |-M| ]
+[ |-Q| ]
 [ |SYN_OPT-R| ]
-[ |-T|\ [*dz*]\ [**+a**\ [*alpha*]]\ [**+s**] ]
+[ |-T|\ [*dv*]\ [**+a**\ [*alpha*]]\ [**+s**] ]
 [ |SYN_OPT-V| ]
 [ |SYN_OPT-f| ]
 [ |SYN_OPT-h| ]
@@ -31,15 +33,16 @@ Synopsis
 Description
 -----------
 
-**grdinfo** reads a 2-D binary grid file and reports metadata and
-various statistics for the (*x*,\ *y*,\ *z*) data in the grid file(s).
-The output information contains the minimum/maximum values for *x*, *y*,
-and *z*, where the min/max of *z* occur, the *x*- and *y*-increments,
-and the number of *x* and *y* nodes, and [optionally] the mean, standard
-deviation, and/or the median, median absolute deviation (MAD) of *z*, and/or
-the mode (Least Median of Squares; LMS), LMS scale of *z*, and number of nodes set
+**grdinfo** reads a 2-D binary grid file and reports metadata and various
+statistics for the data (*v*) and coordinates (*x*,\ *y*) in a grid (or *x*,\ *y*,\ *z* for a 3-D cube).
+The output information may report the minimum/maximum values for *v* and the
+coordinates, where the min/max of *v* occur, the spatial increments,
+and the number of nodes in each dimension, and [optionally] the mean, standard
+deviation, and/or the median, median absolute deviation (MAD) of *v*, and/or
+the mode (Least Median of Squares; LMS), LMS scale of *v*, and number of nodes set
 to NaN. We also report if the grid is pixel- or gridline-registered and
 if it is a Cartesian or Geographic data set (based on metadata in the file).
+With option **-Q** we can also report information for 3-D data cubes.
 
 Required Arguments
 ------------------
@@ -54,16 +57,17 @@ Optional Arguments
 
 **-C**\ [**n**\|\ **t**\]
     Formats the report using tab-separated fields on a single line. The
-    output is *name w e s n z0 z1 dx dy nx ny*\ [ *x0 y0 x1 y1* ] [ *med
+    output is *name w e s n {b t} v0 v1 dx dy {dz} nx ny {nz}*\ [ *x0 y0 {z0} x1 y1 {z1}* ] [ *med
     scale* ] [*mean std rms*] [*n\_nan*] *registration gtype*. The data in brackets are
     output only if the corresponding options **-M**, **-L1**, **-L2**,
-    and **-M** are used, respectively. Use **-Ct** to place file *name*
+    and **-M** are used, respectively, while the data in braces only apply if **-Q** is
+    used with 3-D data cubes. Use **-Ct** to place file *name*
     at the end of the output record or **-Cn** to only output numerical
     columns.  The *registration* is either 0 (gridline) or 1 (pixel),
     while *gtype* is either 0 (Cartesian) or 1 (geographic).
     If the **-I** option is used, the
-    output format is instead *NF w e s n z0 z1*, where *NF* is the total
-    number of grids read and *w e s n* are rounded off (see **-I**).
+    output format is instead *NF w e s n {b t} v0 v1*, where *NF* is the total
+    number of files read and *w e s n {b t}* are rounded off (see **-I**).
 
 .. _-D:
 
@@ -111,15 +115,15 @@ Optional Arguments
 
 **-L**\ [**0**\|\ **1**\|\ **2**\|\ **p**\|\ **a**]
     **-L0**
-        Report range of z after actually scanning the data, not just
+        Report range of *v* after actually scanning the data, not just
         reporting what the header says.
     **-L1**
-        Report median and L1 scale of *z* (L1 scale = 1.4826 \* Median
+        Report median and L1 scale of *v* (L1 scale = 1.4826 \* Median
         Absolute Deviation (MAD)).
     **-L2**
-        Report mean, standard deviation, and root-mean-square (rms) of *z*.
+        Report mean, standard deviation, and root-mean-square (rms) of *v*.
     **-Lp**
-        Report mode (LMS) and LMS scale of *z*.
+        Report mode (LMS) and LMS scale of *v*.
     **-La**
         All of the above.
 
@@ -130,28 +134,35 @@ Optional Arguments
 .. _-M:
 
 **-M**
-    Find and report the location of min/max z-values, and count and
+    Find and report the location of min/max *v*-values, and count and
     report the number of nodes set to NaN, if any.
 
+.. _-Q:
+
+**-Q**
+    All input files must be data 3-D netCDF data cube files [all files are 2-D grids].
+    Not compatible with **-D**, **-E**, **-F**, and **-Ib**.
+
 .. _-R:
 
 .. |Add_-R| replace:: Using the **-R** option
     will select a subsection of the input grid(s). If this subsection
     exceeds the boundaries of the grid, only the common region will be extracted.
+    If **-Q** is used you must also append limits in the *z* dimension.
 .. include:: explain_-R.rst_
 
 .. _-T:
 
-**-T**\ [*dz*]\ [**+a**\ [*alpha*]]\ [**+s**]
-    Determine min and max z-value.  If *dz* is provided then we first round these
-    values off to multiples of *dz*. To exclude the two tails of the distribution
+**-T**\ [*dv*]\ [**+a**\ [*alpha*]]\ [**+s**]
+    Determine min and max data value.  If *dv* is provided then we first round these
+    values off to multiples of *dv*. To exclude the two tails of the distribution
     when determining the min and max you can add **+a** to set the *alpha*
-    value (in percent [2]): We then sort the grid, exclude the data in the
+    value (in percent [2]): We then sort the values, exclude the data in the
     0.5*\ *alpha* and 100 - 0.5*\ *alpha* tails, and revise the min and max.
     To force a symmetrical range about zero, using minus/plus the max
     absolute value of the two extremes, append **+s**. We report the
-    result via the text string **-T**\ *zmin/zmax* or **-T**\ *zmin/zmax/dz*
-    (if *dz* was given) as expected by :doc:`makecpt`.
+    result via the text string **-T**\ *vmin/vmax* or **-T**\ *vmin/vmax/dv*
+    (if *dv* was given) as expected by :doc:`makecpt`.
 
 .. _-V:
 
@@ -181,6 +192,10 @@ Get the grid spacing in earth_relief_10m::
 
     dx=`gmt grdinfo -Cn -o7 @earth_relief_10m`
 
+To learn about the extreme values and coordinates in the 3-D data cube S362ANI_kmps.nc?vs::
+
+    gmt grdinfo -Q -M S362ANI_kmps.nc?vs
+
 See Also
 --------
 

doc/rst/source/grdinterpolate.rst

@@ -12,15 +12,16 @@ Synopsis
 
 .. include:: common_SYN_OPTs.rst_
 
-**gmt grdinterpolate** *3Dgrid* | *grd1 grd2 ...*
+**gmt grdinterpolate** *cube* | *grd1 grd2 ...*
 |-G|\ *outfile*
-|-T|\ [*min/max*\ /]\ *inc*\ [**+i**\|\ **n**] \|\ |-T|\ *file*\|\ *list*
+[ |-D|\ [**+x**\ *xname*][**+y**\ *yname*][**+z**\ *zname*][**+v**\ *vname*][**+s**\ *scale*][**+o**\ *offset*][**+n**\ *invalid*][**+t**\ *title*][**+r**\ *remark*] ]
 [ |-E|\ *table*\|\ *line* ]
 [ |-F|\ **l**\|\ **a**\|\ **c**\|\ **n**\ [**+1**\|\ **2**] ]
 [ |SYN_OPT-R| ]
 [ |-S|\ *x/y*\|\ *pointfile*\ [**+h**\ *header*] ]
+[ |-T|\ [*min/max*\ /]\ *inc*\ [**+i**\|\ **n**] \|\ |-T|\ *file*\|\ *list* ]
 [ |SYN_OPT-V| ]
-[ |-Z|\ **i**\ *levels*\|\ **o** ]
+[ |-Z|\ [*levels*] ]
 [ |SYN_OPT-b| ]
 [ |SYN_OPT-d| ]
 [ |SYN_OPT-e| ]
@@ -42,7 +43,7 @@ Description
 
 **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 *x* and *y* dimensions
+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
@@ -53,28 +54,28 @@ more level-series, or we slice the 3-D cube along an arbitrary vertical slice an
 Required Arguments
 ------------------
 
-*3Dgrid*
-    Name of a 3-D netCDF data cube to be interpolated. Alternatively, with **-Zi**,
+*cube*
+    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.
 
 .. _-G:
 
 **-G**\ *outfile*
     This is the output 3D data cube file.  If **-T** only selects a
     single layer then the data cube collapses to a regular 2-D grid file.
-    If **-Zo** is used then *outfile* must contain a C-format statement
-    for a floating point number.  Also see **-S** for a similar use to
-    write individual level-series tables.
-
-.. _-T:
-
-**-T**\ [*min/max*\ /]\ *inc*\ [**+i**\|\ **n**] \|\ |-T|\ *file*\|\ *list*
-    Make evenly spaced time-steps from *min* to *max* by *inc* [Default uses input times].
-    For details on array creation, see `Generate 1D Array`_.
+    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
 ------------------
 
+.. _-D:
+
+.. include:: explain_-D_cap.rst_
+
 .. _-E:
 
 **-E**\ *table*\|\ *line*\ [,\ *line*,...][**+a**\ *az*][**+g**][**+i**\ *inc*][**+l**\ *length*][**+n**\ *np*][**+o**\ *az*][**+p**][**+r**\ *radius*][**+x**]
@@ -135,21 +136,27 @@ Optional Arguments
     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:
+
+**-T**\ [*min/max*\ /]\ *inc*\ [**+i**\|\ **n**] \|\ |-T|\ *file*\|\ *list*
+    Make evenly spaced time-steps from *min* to *max* by *inc* [Default uses input times].
+    For details on array creation, see `Generate 1D 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.
+
 .. _-V:
 
 .. |Add_-V| unicode:: 0x20 .. just an invisible code
 ..  include:: explain_-V.rst_
 
 .. _-Z:
 
-**-Zi**\ *levels*\|\ **o**
-    Use **-Zi** to obtain the *levels* and then we read the corresponding number of
-    2-D input grids given on the command line [Default is a single 3-D data cube].
-    The *levels* are specified the same way as in **-T**.
-    Use **-Zo** to write the 3-D data cube as a series of 2-D grids instead.  If used,
-    then the *outgrid* name given by **-G** must contain a C-language format statement
-    for a floating point number (for instance, one can try layer_%6.6f.grd) which will contain the level
-    for each grid [Default is a 3-D data cube, unless only one layer is implied by **-T**].
+**-Z**\ [*levels*]
+    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)
@@ -194,8 +201,8 @@ File Order
 ----------
 
 If you provide a series of 2-D files and thus separately assigning the
-level via **-Zi**, then you must make sure that the order the grids are given
-on the command line matches the levels you provide via **-Zi**.  Unless your
+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).
 
@@ -234,7 +241,7 @@ To extract a single, new 2-D layer from the 3-D cube implied by the individual g
 layers_*.nc, with individual layer values given via z.txt, for level 3400
 using a linear spline, try::
 
-    gmt grdinterpolate layers_*.nc -Ziz.txt -T3400 -Fl -Gtemp_3400.nc
+    gmt grdinterpolate layers_*.nc -Zz.txt -T3400 -Fl -Gtemp_3400.nc
 
 To resample the the temperature.nc 3-D cube for all levels from
 1500 to 2500 in steps of 50, using an Akima spline, try::
@@ -243,13 +250,13 @@ To resample the the temperature.nc 3-D cube for all levels from
 
 The same, but this time write individual 2-D grids per layer, try::
 
-    gmt grdinterpolate temperature.nc -T1500/2500/50 -Gtemperature_%4.0f.nc -Fa -Zo
+    gmt grdinterpolate temperature.nc -T1500/2500/50 -Gtemperature_%4.0f.nc -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 -Zidates.txt -S115W/33N+h"Some like it hot" > record.txt
+    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 S362ANI_kmps.nc with seismic velocities that goes
 through the Hawaii hotspot, selecting cube vs (Isotropic Shear Velocity) and letting the