grdgradient

Compute directional gradients from a grid

Synopsis

gmt grdgradient ingrid -Goutgrid [ -Aazim[/azim2] ] [ -D[a][c][n][o] ] [ -E[m|s|p]azim/elev[+aambient][+ddiffuse][+pspecular][+sshine] ] [ -N[e|t][amp][+aambient][+ssigma][+ooffset] ] [ -Qc|r|R[+ffile]] [ -Rregion ] [ -Sslopefile ] [ -V[level] ] [ -fflags ] [ -nflags ] [ --PAR=value ]

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

Description

grdgradient may be used to compute the directional derivative in a given direction (-A), or to find the direction (-D) [and the magnitude (-S)] of the vector gradient of the data.

Estimated values in the first/last row/column of output depend on boundary conditions (see -n).

Required Arguments

ingrid[=ID|?varname][+bband][+ddivisor][+ninvalid][+ooffset][+sscale]

Optionally, append =ID for reading a specific file format [Default is =nf] or ?varname for a specific netCDF variable [Default is the first 2-D grid found by GMT]. The following modifiers are supported:

  • +b - Select a band [Default is 0].

  • +d - Divide data values by the given divisor [Default is 1].

  • +n - Replace data values matching invalid with NaN.

  • +o - Offset data values by the given offset [Default is 0].

  • +s - Scale data values by the given scale [Default is 1].

Note: Any offset is added after any scaling.

-Goutgrid[=ID][+ddivisor][+ninvalid][+ooffset|a][+sscale|a][:driver[dataType][+coptions]]

Optionally, append =ID for writing a specific file format. The following modifiers are supported:

  • +d - Divide data values by given divisor [Default is 1].

  • +n - Replace data values matching invalid with a NaN.

  • +o - Offset data values by the given offset, or append a for automatic range offset to preserve precision for integer grids [Default is 0].

  • +s - Scale data values by the given scale, or append a for automatic scaling to preserve precision for integer grids [Default is 1].

Note: Any offset is added before any scaling. +sa also sets +oa (unless overridden). To write specific formats via GDAL, use =gd and supply driver (and optionally dataType) and/or one or more concatenated GDAL -co options using +c. See the “Writing grids and images” cookbook section for more details.

Optional Arguments

-Aazim[/azim2]

Azimuthal direction for a directional derivative; azim is the angle in the x-y plane measured in degrees positive clockwise from north (the +y direction) toward east (the +x direction). The negative of the directional derivative, \(-(\frac{dz}{dx}\sin(a) + \frac{dz}{dy}\cos(a))\) , is found where \(a\) is the azimuth; negation yields positive values when the slope of \(z(x,y)\) is downhill in the \(a\) direction, the correct sense for shading the illumination of an image (see grdimage and grdview) by a light source above the x-y plane shining from the \(a\) direction. Optionally, supply two azimuths, -Aazim/azim2, in which case the gradients in each of these directions are calculated and the one larger in magnitude is retained; this is useful for illuminating data with two directions of lineated structures, e.g., -A0/270 illuminates from the north (top) and west (left). Finally, if azim is a file it must be a grid of the same domain, spacing and registration as ingrid and we will update the azimuth at each output node when computing the directional derivatives.

-D[a][c][n][o]

Instead of finding the gradient in a given direction (-A), find the direction of the positive (up-slope) gradient of the data. The following directives are supported:

  • a: Find the aspect (i.e., the down-slope direction).

  • c: Use the conventional Cartesian angles measured counterclockwise from the positive x (east direction) [Default measures directions clockwise from north, as azim in -A].

  • n: Add 90 degrees to all angles (e.g., to give local strikes of the surface).

  • o: Report orientations (0-180) rather than directions (0-360).

-E[m|s|p]azim/elev[+aambient][+ddiffuse][+pspecular][+sshine]

Compute Lambertian radiance appropriate to use with grdimage and grdview. The Lambertian Reflection assumes an ideal surface that reflects all the light that strikes it and the surface appears equally bright from all viewing directions. Here, azim and elev are the azimuth and elevation of the light vector. Optionally, supply ambient [0.55], diffuse [0.6], specular [0.4], or shine [10], which are parameters that control the reflectance properties of the surface. Default values are given in the brackets. Use -Es for a simpler Lambertian algorithm. Note that with this form you only have to provide azimuth and elevation. Alternatively, use -Ep for the Peucker piecewise linear approximation (simpler but faster algorithm; in this case the azim and elev are hardwired to 315 and 45 degrees. This means that even if you provide other values they will be ignored.)

-N[e|t][amp][+aambient][+ssigma][+ooffset]

Normalization [Default is no normalization.] The actual gradients \(g\) are offset and scaled to produce normalized gradients \(g_n\) with a maximum output magnitude of amp. If amp is not given, default amp = 1. If offset is not given, it is set to the average of \(g\). The following forms are supported, where \(o\) is the offset and \(a\) is the amp:

  • -N - Normalize using \(g_n = a(\frac{g - o}{max(|g - o|)})\)

  • -Ne - Normalize using a cumulative Laplace distribution yielding: \(g_n = a(1 - \exp{(\sqrt{2}\frac{g - o}{\sigma}))}\), where \(\sigma\) is estimated using the L1 norm of \((g - o)\) if it is not given.

  • -Nt - Normalize using a cumulative Cauchy distribution yielding: \(g_n = \frac{2a}{\pi}(\tan^{-1}(\frac{g - o}{\sigma}))\) where \(\sigma\) is estimated using the L2 norm of \((g -o)\) if it is not given.

To use \(o\) and/or \(\sigma\) from a previous calculation, leave out the argument to the modifier(s) and see -Q for usage. As a final option, you may add +aambient to add ambient to all nodes after gradient calculations are completed.

-Qc|r|R[+ffile]

Controls how normalization via -N is carried out. When multiple grids should be normalized the same way (i.e., with the same offset and/or sigma), we must pass these values via -N. However, this is inconvenient if we compute these values from a grid. Use -Qc to save the results of offset and sigma to a statistics file; if grid output is not needed for this run then do not specify -G. For subsequent runs, just use -Qr to read these values. Using -QR will read then delete the statistics file. See Tiles for more information. Optionally, append +ffile to write/read the statistics to/from the specified file [Default is grdgradient.stat in the current TMP directory].

-Rxmin/xmax/ymin/ymax[+r][+uunit]

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

-Sslopefile

Name of output grid file with scalar magnitudes of gradient vectors. Requires -D but makes -G optional.

-V[level]

Select verbosity level [w]. (See full description) (See technical reference).

-fflags

Geographic grids (dimensions of longitude, latitude) will be converted to meters via a “Flat Earth” approximation using the current ellipsoid parameters.

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

Select interpolation mode for grids.

-^ 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.

--PAR=value

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

Grid Distance Units

If the grid does not have meter as the horizontal unit, append +uunit to the input file name to convert from the specified unit to meter. If your grid is geographic, convert distances to meters by supplying -fflags instead.

Hints

  • If you don’t know what -N options to use to make an intensity file for grdimage or grdview, a good first try is -Ne0.6.

  • Usually 255 shades are more than enough for visualization purposes. You can save 75% disk space by appending =nb/a to the output filename outgrid.

  • If you want to make several illuminated maps of subregions of a large data set, and you need the illumination effects to be consistent across all the maps, use the -N option and supply the same value of sigma and offset to grdgradient for each map. A good guess is offset = 0 and sigma found by grdinfo -L2 or -L1 applied to an unnormalized gradient grd.

  • If you simply need the x- or y-derivatives of the grid, use grdmath.

Tiles

For very large datasets (or very large plots) you may need to break the job into multiple tiles. It is then important that the normalization of the intensities are handled the same way for each tile. By default, offset and sigma are recalculated for each tile. Hence, different tiles of the same large grid will compute different offset and sigma values. Thus, the intensity for the same directional slope will be different across the final map. This inconsistency can lead to visible changes in image appearance across tile seams. The way to ensure compatible results is to specify the same offset and sigma via the modifiers to -N. However, if these need to be estimated from the large grid then the -Q option can help: Run grdgradient on the full grid (or as large portion of the grid that your computer can handle) and specify -Qc to create a statistics file with the resulting offset and sigma. Then, for each of your grid tile calculations, give +o and/or +s without arguments to -N and specify -Qr. This option will read the values from the hidden statistics file and use them in the normalization. If you use -QR for the final tile then the statistics file is removed after use.

Ambient

The ambient light offset is used to darken or brighten all intensities. This modifier is typically used to darken an entire image by subtracting a constant from all the intensities. E.g., if you use +a-0.5 then you subtract 0.5 from all intensities, making them more negative and hence darken the image.

Special Effects

Users who wishes to highlight just an area of a grid are encouraged to consider creative uses of the vast number of operators in grdmath and the grid clipping schemes in grdmask and grdclip. Many separate intensity grids can be combined into one via grid algebra and boolean operators in grdmath. The sky is the limit.

Examples

Note: Below are some examples of valid syntax for this module. The examples that use remote files (file names starting with @) can be cut and pasted into your terminal for testing. Other commands requiring input files are just dummy examples of the types of uses that are common but cannot be run verbatim as written.

To make a file for illuminating a portion of the data in the remote file @earth_relief_05 using exponentially normalized gradients in the range [-0.6,0.6] imitating light sources in the north and west directions:

gmt grdgradient @earth_relief_05m -R0/20/0/20 -A0/270 -Ggradients.nc -Ne0.6 -V

To find the azimuth orientations of seafloor fabric in the file topo.nc:

gmt grdgradient topo.nc -Dno -Gazimuths.nc -V

To determine the offset and sigma suitable for normalizing the intensities from topo.nc, do

gmt grdgradient topo.nc -A30 -Nt0.6 -Qc -V

Without -G, only the hidden statistics file is created and no output grid is written.

To use the previously determined offset and sigma to normalize the intensities in tile_3.nc, do

gmt grdgradient tile_3.nc -A30 -Nt0.6+o+s -Qr -V -Gtile_3_int.nc

References

Horn, B.K.P., Hill-Shading and the Reflectance Map, Proceedings of the IEEE, Vol. 69, No. 1, January 1981, pp. 14-47. (https://www.researchgate.net/publication/2996084_Hill_shading_and_the_reflectance_map)

See Also

gmt, gmt.conf grdhisteq, grdinfo, grdmath, grdimage, grdview, grdvector