gpsgridder

Interpolate GPS velocities using Green’s functions for elastic deformation

Synopsis

gmt gpsgridder [ table ] -Goutgrid [ -C[[n|r|v]value[%]][+c][+ffile][+i][+n] ] [ -E[misfitfile][+rreportfile] ] [ -F[d|f]fudge ] [ -Iincrement ] [ -L ] [ -Nnodefile ] [ -Rregion ] [ -Snu ] [ -Tmaskgrid ] [ -V[level] ] [ -W[+s|w] ] [ -bbinary ] [ -dnodata[+ccol] ] [ -eregexp ] [ -fflags ] [ -hheaders ] [ -oflags ] [ -qiflags ] [ -x[[-]n] ] [ -:[i|o] ] [ --PAR=value ]

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

Description

gpsgridder grids 2-D vector data such as GPS velocities by using a coupled model based on 2-D elasticity. The degree of coupling can be tuned by adjusting the effective Poisson’s ratio, \(\nu\). The solution field can be tuned to extremes such as incompressible (1), typical elastic (0.5) or even an unphysical value of -1 that basically removes the elastic coupling of vector interpolation. Smoothing is offered via the optional elimination of small eigenvalues. The solutions for the two component grids are evaluated as

\[\begin{split}u(\mathbf{x}) = \sum_{j=1}^{n} \alpha_j q(\mathbf{x}, \mathbf{x}_j) + \beta_j w(\mathbf{x}, \mathbf{x}_j)\\ v(\mathbf{x}) = \sum_{j=1}^{n} \alpha_j w(\mathbf{x}, \mathbf{x}_j) + \beta_j p(\mathbf{x}, \mathbf{x}_j)\end{split}\]

where the three 2-D elastic coupled Green’s functions are given by

\[\begin{split}q(\mathbf{a}, \mathbf{b}) = (3 - \nu)\log r + (1 + \nu) \frac{y^2}{r^2}\\ p(\mathbf{a}, \mathbf{b}) = (3 - \nu)\log r + (1 + \nu) \frac{x^2}{r^2}\\ w(\mathbf{a}, \mathbf{b}) = -(1 + \nu) \frac{xy}{r^2}\end{split}\]

Here, r is the radial distance between points a and b and x and y are the components of that distance. The body forces \(\alpha_j\) and \(\beta_j\) are obtained by evaluating the solution at the data locations and inverting the square linear system that results; see Sandwell and Wessel [2016] and Haines et al. [2015] for details.

Required Arguments

table

table with GPS strain rates at discrete locations. We expect the input format to be x y u v [ du dv ] (see -W to specify data uncertainties or weights). If lon lat is given you must supply -fg and we will use a flat Earth approximation in the calculation of distances.

-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

-C[[n|r|v]value[%]][+c][+ffile][+i][+n]
Find an approximate surface fit: Solve the linear system for the

spline coefficients by SVD and eliminate the contribution from smaller eigenvalues [Default uses Gauss-Jordan elimination to solve the linear system and fit the data exactly (unless -W is used)]. Append a directive and value to determine which eigenvalues to keep:

  • n - Retain only the value numbers largest eigenvalues [all]. Optionally, append % to indicate value is given in percentage.

  • r - Retain those eigenvalues whose ratio to the largest eigenvalue is less than value [Default, with value = 0].

  • v - Retain the eigenvalues needed to ensure the model prediction variance fraction is at least value. Optionally, append % to indicate value is given in percentage.

Several optional modifiers are available:

  • +c - Produce the cumulative sum of these contributions, two (u and v) grids per eigenvalue (2-D only).

  • +f - Append file to save the eigenvalues to the specified file for further analysis.

  • +i - Produce the incremental sum of these contributions, two (u and v) grids per eigenvalue (2-D only).

  • +n - If given then +ffile is required and execution will stop after saving the eigenvalues, i.e., no surface output is produced.

Notes: (1) Modifiers ++c and +i require a file name with a suitable extension to be given via -G (we automatically insert “_cum_###” or “_inc_###” before the extension, using a fixed integer format for the eigenvalue number, starting at 0). (2) Use both modifiers to write both types of intermediate grids.

-E[misfitfile][+rreportfile]

Evaluate the spline exactly at the input data locations and report statistics of the misfit (mean, standard deviation, and rms) for u and v separately and combined. Optionally, append a filename and we will write the data table, augmented by two extra columns after each of the u and v columns holding the spline estimates and misfits. If -W is given we also add two more columns with \(\chi_u^2\) and \(\chi_v^2\) values. Alternatively, if -C is used and history is computed (via one or more of modifiers +c and +i), then we will instead write a table with eigenvalue number, eigenvalue, percent of model variance explained, and overall rms, rms_u, and rms_v misfits. If -W is used we also append \(\chi^2\), \(\chi_u^2\), and \(\chi_v^2\).

Optionally append this modifier:

  • +r - Write misfit and variance evaluation statistics to reportfile. Output order is Data Model Explained(%) N Mean Std.dev RMS. If -W is used we add \(\chi^2\) as a final 8th column. Note: If +r is not used then -Vi will report the information via verbose messages.

-F[d|f]fudge

The Green’s functions are proportional to \(r^{-2}\) and \(\log(r)\) and thus blow up for r == 0. To prevent that we offer two fudging schemes: -Fddel_radius lets you add a constant offset to all radii and must be specified in the user units. Alternatively, use -Fffactor which will compute del_radius from the product of the shortest inter-point distance and factor [0.01].

-Ix_inc[+e|n][/y_inc[+e|n]]

Set the grid spacing as x_inc [and optionally y_inc].

Geographical (degrees) coordinates: Optionally, append an increment unit. Choose among:

  • d - Indicate arc degrees

  • m - Indicate arc minutes

  • s - Indicate arc seconds

If one of e (meter), f (foot), k (km), M (mile), n (nautical mile) or u (US survey foot), the increment will be converted to the equivalent degrees longitude at the middle latitude of the region (the conversion depends on PROJ_ELLIPSOID). If y_inc is not given or given but set to 0 it will be reset equal to x_inc; otherwise it will be converted to degrees latitude.

All coordinates: The following modifiers are supported:

  • +e - Slightly adjust the max x (east) or y (north) to fit exactly the given increment if needed [Default is to slightly adjust the increment to fit the given domain].

  • +n - Define the number of nodes rather than the increment, in which case the increment is recalculated from the number of nodes, the registration (see GMT File Formats), and the domain. Note: If -Rgrdfile is used then the grid spacing and the registration have already been initialized; use -I and -R to override these values.

-L

Leave trend alone. Do not remove a planer (2-D) trend from the data before fitting the spline. [Default removes least squares plane, fits normalized residuals, and restores plane].

-Nnodefile

ASCII file with coordinates of desired output locations x in the first column(s). The resulting w values are appended to each record and written to the file given in -G [or standard output if not specified]; see -bo for binary output instead. This option eliminates the need to specify options -R, -I, and -r.

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

Specify the region of interest. (See full description) (See technical reference).

-Snu

Specify Poisson’s ratio to use for this 2-D elastic sheet [0.5]. Note: 1.0 is incompressible in a 2-D formulation while -1 removes all coupling between the two directions.

-Tmaskgrid

Only evaluate the solutions at the nodes in the maskgrid that are not set to NaN. This option eliminates the need to specify options -R, -I (and -r).

-W[+s|w]

One-sigma data uncertainties for u and v are provided in the last two columns. We then compute least squares weights that are inversely proportional to the square of the uncertainties [Default, or +s]. Instead, append +w if weights are given instead of uncertainties, in which case we just use the weights as provided (no squaring). This results in a weighted least squares fit. Note that -W only has an effect if -C is used [Default uses no weights or uncertainties].

-V[level]

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

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

-fg

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

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

Skip or produce header record(s). Not used with binary data.

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

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

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

-r[g|p] (more …)

Set node registration [gridline].

-:[i|o] (more …)

Swap 1st and 2nd column on input and/or 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.

--PAR=value

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

Units

For map distance unit, append unit d for arc degree, m for arc minute, and s for arc second, or e for meter [Default unless stated otherwise], f for foot, k for km, M for statute mile, n for nautical mile, and u for US survey foot. By default we compute such distances using a spherical approximation with great circles (-jg) using the authalic radius (see PROJ_MEAN_RADIUS). You can use -jf to perform “Flat Earth” calculations (quicker but less accurate) or -je to perform exact geodesic calculations (slower but more accurate; see PROJ_GEODESIC for method used).

Notes on SVD solution

It may be difficult to know how many eigenvalues are needed for a suitable approximate fit. The -C modifiers allow you to explore this further by creating solutions for all cutoff selections and estimate model variance and data misfit as a function of how many eigenvalues are used. The large set of such solutions can be animated so it is easier to explore the changes between solutions and to make a good selection for the -C directive values. See the animations for one or more examples of this exploration.

Examples

To compute the u and v strain rate grids from the GPS data set gps.txt, containing x y u v du dv, on a 2x2 arc minute grid for California, and just using about 25% of the largest eigenvalues, try:

gmt gpsgridder gps.txt -R-125/-114/31/41 -I2m -fg -W -r -Cn25% -Ggps_strain_%s.nc -V

Deprecations

  • 6.3.0: Use +n instead of negative value for -C to set dry-run. #5725

References

Haines, A. J. et al., 2015, Enhanced Surface Imaging of Crustal Deformation, SpringerBriefs in Earth Sciences, https://doi.org/10.1007/978-3-319-21578-5_2.

Sandwell, D. T. and P. Wessel, 2016, Interpolation of 2-D Vector Data Using Constraints from Elasticity, Geophys. Res. Lett., 43, 10,703-10,709, https://doi.org/10.1002/2016GL070340

See Also

gmt, greenspline nearneighbor, surface