pygmt.grdtrack

pygmt.grdtrack(grid, points=None, newcolname=None, outfile=None, *, resample=None, crossprofile=None, dfile=None, profile=None, critical=None, region=None, no_skip=None, stack=None, radius=None, verbose=None, z_only=None, aspatial=None, binary=None, nodata=None, find=None, coltypes=None, gap=None, header=None, incols=None, distcalc=None, interpolation=None, outcols=None, skiprows=None, wrap=None, **kwargs)

Sample grids at specified (x,y) locations.

Reads one or more grid files and a table (from file or an array input; but see profile for exception) with (x,y) [or (lon,lat)] positions in the first two columns (more columns may be present). It interpolates the grid(s) at the positions in the table and writes out the table with the interpolated values added as (one or more) new columns. Alternatively (crossprofile), the input is considered to be line-segments and we create orthogonal cross-profiles at each data point or with an equidistant separation and sample the grid(s) along these profiles. A bicubic [Default], bilinear, B-spline or nearest-neighbor interpolation is used, requiring boundary conditions at the limits of the region (see interpolation; Default uses “natural” conditions (second partial derivative normal to edge is zero) unless the grid is automatically recognized as periodic.)

Full option list at https://docs.generic-mapping-tools.org/latest/grdtrack.html

Aliases:

• A = resample

• C = crossprofile

• D = dfile

• E = profile

• F = critical

• N = no_skip

• R = region

• S = stack

• T = radius

• V = verbose

• Z = z_only

• a = aspatial

• b = binary

• d = nodata

• e = find

• f = coltypes

• g = gap

• h = header

• i = incols

• j = distcalc

• n = interpolation

• o = outcols

• s = skiprows

• w = wrap

Parameters

• grid (xarray.DataArray or str) – Gridded array from which to sample values from, or a filename (netcdf format).

• points (str or numpy.ndarray or pandas.DataFrame or xarray.Dataset or geopandas.GeoDataFrame) – Pass in either a file name to an ASCII data table, a 2D numpy.ndarray, a pandas.DataFrame, an xarray.Dataset made up of 1D xarray.DataArray data variables, or a geopandas.GeoDataFrame containing the tabular data.

• newcolname (str) – Required if points is a pandas.DataFrame. The name for the new column in the track pandas.DataFrame table where the sampled values will be placed.

• outfile (str) – The file name for the output ASCII file.

• resample (str) – f|p|m|r|R[+l] For track resampling (if crossprofile or profile are set) we can select how this is to be performed. Append f to keep original points, but add intermediate points if needed [Default], m as f, but first follow meridian (along y) then parallel (along x), p as f, but first follow parallel (along y) then meridian (along x), r to resample at equidistant locations; input points are not necessarily included in the output, and R as r, but adjust given spacing to fit the track length exactly. Finally, append +l if geographic distances should be measured along rhumb lines (loxodromes) instead of great circles. Ignored unless crossprofile is used.

• crossprofile (str) – length/ds[/spacing][+a|+v][l|r]. Use input line segments to create an equidistant and (optionally) equally-spaced set of crossing profiles along which we sample the grid(s) [Default simply samples the grid(s) at the input locations]. Specify two length scales that control how the sampling is done: length sets the full length of each cross-profile, while ds is the sampling spacing along each cross-profile. Optionally, append /spacing for an equidistant spacing between cross-profiles [Default erects cross-profiles at the input coordinates]; see resample for how resampling the input track is controlled. By default, all cross-profiles have the same direction (left to right as we look in the direction of the input line segment). Append +a to alternate the direction of cross-profiles, or v to enforce either a “west-to-east” or “south-to-north” view. By default the entire profiles are output. Choose to only output the left or right halves of the profiles by appending +l or +r, respectively. Append suitable units to length; it sets the unit used for ds [and spacing] (See Units). The default unit for geographic grids is meter while Cartesian grids implies the user unit. The output columns will be lon, lat, dist, azimuth, z1, z2, …, zn (The zi are the sampled values for each of the n grids).

• dfile (str) – In concert with crossprofile we can save the (possibly resampled) original lines to dfile [Default only saves the cross-profiles]. The columns will be lon, lat, dist, azimuth, z1, z2, … (sampled value for each grid).

• profile (str) – line[,line,…][+aaz][+c][+d][+g][+iinc][+llength][+nnp][+oaz][+rradius]. Instead of reading input track coordinates, specify profiles via coordinates and modifiers. The format of each line is start/stop, where start or stop are either lon/lat (x/y for Cartesian data) or a 2-character XY key that uses the text-style justification format to specify a point on the map as [LCR][BMT]. Each line will be a separate segment unless +c is used which will connect segments with shared joints into a single segment. In addition to line coordinates, you can use Z-, Z+ to mean the global minimum and maximum locations in the grid (only available if a single grid is given via outfile). You may append +iinc to set the sampling interval; if not given then we default to half the minimum grid interval. 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. 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. Append +d to output the along-track distances after the coordinates. Note: No track file will be read. 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 distcalc to control distance calculation mode [Default is Great Circle]. Note: If crossprofile is set and spacing is given then that sampling scheme overrules any modifier set in profile.

• critical (str) – [+b][+n][+r][+zz0]. Find critical points along each cross-profile as a function of along-track distance. Requires crossprofile and a single input grid (z). We examine each cross-profile generated and report (dist, lonc, latc, distc, azimuthc, zc) at the center peak of maximum z value, (lonl, latl, distl) and (lonr, latr, distr) at the first and last non-NaN point whose z-value exceeds z0, respectively, and the width based on the two extreme points found. Here, dist is the distance along the original input points and the other 12 output columns are a function of that distance. When searching for the center peak and the extreme first and last values that exceed the threshold we assume the profile is positive up. If we instead are looking for a trough then you must use +n to temporarily flip the profile to positive. The threshold z0 value is always given as >= 0; use +z to change it [Default is 0]. Alternatively, use +b to determine the balance point and standard deviation of the profile; this is the weighted mean and weighted standard deviation of the distances, with z acting as the weight. Finally, use +r to obtain the weighted rms about the cross-track center (distc == 0). Note: We round the exact results to the nearest distance nodes along the cross-profiles. We write 13 output columns per track: dist, lonc, latc, distc, azimuthc, zc, lonl, latl, distl, lonr, latr, distr, width.

• region (str or list) – xmin/xmax/ymin/ymax[+r][+uunit]. Specify the region of interest.

• no_skip (bool) – Do not skip points that fall outside the domain of the grid(s) [Default only output points within grid domain].

• stack (str or list) –

  method/modifiers. In conjunction with crossprofile, compute a single stacked profile from all profiles across each segment. Choose how stacking should be computed [Default method is a]:

  • a = mean (average)

  • m = median

  • p = mode (maximum likelihood)

  • l = lower

  • L = lower but only consider positive values

  • u = upper

  • U = upper but only consider negative values.

  The modifiers control the output; choose one or more among these choices:

  • +a : Append stacked values to all cross-profiles.

  • +d : Append stack deviations to all cross-profiles.

  • +r : Append data residuals (data - stack) to all cross-profiles.

  • +s[file] : Save stacked profile to file [Default filename is grdtrack_stacked_profile.txt].

  • +cfact : Compute envelope on stacked profile as ±fact *deviation [Default fact value is 2].

  Notes:

  1. Deviations depend on method and are st.dev (a), L1 scale, i.e., 1.4826 * median absolute deviation (MAD) (for m and p), or half-range (upper-lower)/2.

  2. The stacked profile file contains a leading column plus groups of 4-6 columns, with one group for each sampled grid. The leading column holds cross distance, while the first four columns in a group hold stacked value, deviation, min value, and max value, respectively. If method is one of a|m|p then we also write the lower and upper confidence bounds (see +c). When one or more of +a, +d, and +r are used then we also append the stacking results to the end of each row, for all cross-profiles. The order is always stacked value (+a), followed by deviations (+d) and finally residuals (+r). When more than one grid is sampled this sequence of 1-3 columns is repeated for each grid.

• radius (bool or int or float or str) – [radius][+e|p]. To be used with normal grid sampling, and limited to a single, non-IMG grid. If the nearest node to the input point is NaN, search outwards until we find the nearest non-NaN node and report that value instead. Optionally specify a search radius which limits the consideration to points within this distance from the input point. To report the location of the nearest node and its distance from the input point, append +e. The default unit for geographic grid distances is spherical degrees. Use radius to change the unit and give radius = 0 if you do not want to limit the radius search. To instead replace the input point with the coordinates of the nearest node, append +p.

• verbose (bool or str) –

  Select verbosity level [Default is w], which modulates the messages written to stderr. Choose among 7 levels of verbosity:

  • q - Quiet, not even fatal error messages are produced

  • e - Error messages only

  • w - Warnings [Default]

  • t - Timings (report runtimes for time-intensive algorithms);

  • i - Informational messages (same as verbose=True)

  • c - Compatibility warnings

  • d - Debugging messages

• z_only (bool) – Only write out the sampled z-values [Default writes all columns].

• aspatial (bool or str) – [col=]name[,…]. Control how aspatial data are handled during input and output. Full documentation is at https://docs.generic-mapping-tools.org/latest/gmt.html#aspatial-full.

• binary (bool or str) –

  i|o[ncols][type][w][+l|b]. Select native binary input (using binary="i") or output (using binary="o"), where ncols is the number of data columns of type, which must be one of:

  • c - int8_t (1-byte signed char)

  • u - uint8_t (1-byte unsigned char)

  • h - int16_t (2-byte signed int)

  • H - uint16_t (2-byte unsigned int)

  • i - int32_t (4-byte signed int)

  • I - uint32_t (4-byte unsigned int)

  • l - int64_t (8-byte signed int)

  • L - uint64_t (8-byte unsigned int)

  • f - 4-byte single-precision float

  • d - 8-byte double-precision float

  • x - use to skip ncols anywhere in the record

  For records with mixed types, append additional comma-separated combinations of ncols type (no space). The following modifiers are supported:

  • w after any item to force byte-swapping.

  • +l|b to indicate that the entire data file should be read as little- or big-endian, respectively.

  Full documentation is at https://docs.generic-mapping-tools.org/latest/gmt.html#bi-full.

• nodata (str) – i|onodata. Substitute specific values with NaN (for tabular data). For example, d="-9999" will replace all values equal to -9999 with NaN during input and all NaN values with -9999 during output. Prepend i to the nodata value for input columns only. Prepend o to the nodata value for output columns only.

• find (str) – [~]“pattern” | [~]/regexp/[i]. Only pass records that match the given pattern or regular expressions [Default processes all records]. Prepend ~ to the pattern or regexp to instead only pass data expressions that do not match the pattern. Append i for case insensitive matching. This does not apply to headers or segment headers.

• coltypes (str) – [i|o]colinfo. Specify data types of input and/or output columns (time or geographical data). Full documentation is at https://docs.generic-mapping-tools.org/latest/gmt.html#f-full.

• gap (str or list) –

  x|y|z|d|X|Y|Dgap[u][+a][+ccol][+n|p]. Examine the spacing between consecutive data points in order to impose breaks in the line. To specify multiple criteria, provide a list with each item containing a string describing one set of criteria.

  • x|X - define a gap when there is a large enough change in the x coordinates (upper case to use projected coordinates).

  • y|Y - define a gap when there is a large enough change in the y coordinates (upper case to use projected coordinates).

  • d|D - define a gap when there is a large enough distance between coordinates (upper case to use projected coordinates).

  • z - define a gap when there is a large enough change in the z data. Use +ccol to change the z data column [Default col is 2 (i.e., 3rd column)].

  A unit u may be appended to the specified gap:

  • For geographic data (x|y|d), the unit may be arc d(egree), m(inute), and s(econd), or (m)e(ter), f(eet), k(ilometer), M(iles), or n(autical miles) [Default is (m)e(ter)].

  • For projected data (X|Y|D), the unit may be i(nch), c(entimeter), or p(oint).

  Append modifier +a to specify that all the criteria must be met [default imposes breaks if any one criterion is met].

  One of the following modifiers can be appended:

  • +n - specify that the previous value minus the current column value must exceed gap for a break to be imposed.

  • +p - specify that the current value minus the previous value must exceed gap for a break to be imposed.

• header (str) –

  [i|o][n][+c][+d][+msegheader][+rremark][+ttitle]. Specify that input and/or output file(s) have n header records [Default is 0]. Prepend i if only the primary input should have header records. Prepend o to control the writing of header records, with the following modifiers supported:

  • +d to remove existing header records.

  • +c to add a header comment with column names to the output [Default is no column names].

  • +m to add a segment header segheader to the output after the header block [Default is no segment header].

  • +r to add a remark comment to the output [Default is no comment]. The remark string may contain \n to indicate line-breaks.

  • +t to add a title comment to the output [Default is no title]. The title string may contain \n to indicate line-breaks.

  Blank lines and lines starting with # are always skipped.

• incols (str or 1d array) –

  Specify data columns for primary input in arbitrary order. Columns can be repeated and columns not listed will be skipped [Default reads all columns in order, starting with the first (i.e., column 0)].

  • For 1d array: specify individual columns in input order (e.g., incols=[1,0] for the 2nd column followed by the 1st column).

  • For str: specify individual columns or column ranges in the format start[:inc]:stop, where inc defaults to 1 if not specified, with columns and/or column ranges separated by commas (e.g., incols="0:2,4+l" to input the first three columns followed by the log-transformed 5th column). To read from a given column until the end of the record, leave off stop when specifying the column range. To read trailing text, add the column t. Append the word number to t to ingest only a single word from the trailing text. Instead of specifying columns, use incols="n" to simply read numerical input and skip trailing text. Optionally, append one of the following modifiers to any column or column range to transform the input columns:

    • +l to take the log10 of the input values.

    • +d to divide the input values by the factor divisor [Default is 1].

    • +s to multiple the input values by the factor scale [Default is 1].

    • +o to add the given offset to the input values [Default is 0].

• distcalc (str) –

  e|f|g. Determine how spherical distances are calculated.

  • e - Ellipsoidal (or geodesic) mode

  • f - Flat Earth mode

  • g - Great circle distance [Default]

  All spherical distance calculations depend on the current ellipsoid (PROJ_ELLIPSOID), the definition of the mean radius (PROJ_MEAN_RADIUS), and the specification of latitude type (PROJ_AUX_LATITUDE). Geodesic distance calculations is also controlled by method (PROJ_GEODESIC).

• interpolation (str) –

  [b|c|l|n][+a][+bBC][+c][+tthreshold]. Select interpolation mode for grids. You can select the type of spline used:

  • b for B-spline

  • c for bicubic [Default]

  • l for bilinear

  • n for nearest-neighbor

• outcols (str or 1d array) –

  cols[,…][,t[word]]. Specify data columns for primary output in arbitrary order. Columns can be repeated and columns not listed will be skipped [Default writes all columns in order, starting with the first (i.e., column 0)].

  • For 1d array: specify individual columns in output order (e.g., outcols=[1,0] for the 2nd column followed by the 1st column).

  • For str: specify individual columns or column ranges in the format start[:inc]:stop, where inc defaults to 1 if not specified, with columns and/or column ranges separated by commas (e.g., outcols="0:2,4" to output the first three columns followed by the 5th column). To write from a given column until the end of the record, leave off stop when specifying the column range. To write trailing text, add the column t. Append the word number to t to write only a single word from the trailing text. Instead of specifying columns, use outcols="n" to simply read numerical input and skip trailing text. Note: if incols is also used then the columns given to outcols correspond to the order after the incols selection has taken place.

• skiprows (bool or str) –

  [cols][+a][+r]. Suppress output for records whose z-value equals NaN [Default outputs all records]. Optionally, supply a comma-separated list of all columns or column ranges to consider for this NaN test [Default only considers the third data column (i.e., cols = 2)]. Column ranges must be given in the format start[:inc]:stop, where inc defaults to 1 if not specified. The following modifiers are supported:

  • +r to reverse the suppression, i.e., only output the records whose z-value equals NaN.

  • +a to suppress the output of the record if just one or more of the columns equal NaN [Default skips record only if values in all specified cols equal NaN].

• wrap (str) –

  y|a|w|d|h|m|s|cperiod[/phase][+ccol]. Convert the input x-coordinate to a cyclical coordinate, or a different column if selected via +ccol. The following cyclical coordinate transformations are supported:

  • y - yearly cycle (normalized)

  • a - annual cycle (monthly)

  • w - weekly cycle (day)

  • d - daily cycle (hour)

  • h - hourly cycle (minute)

  • m - minute cycle (second)

  • s - second cycle (second)

  • c - custom cycle (normalized)

  Full documentation is at https://docs.generic-mapping-tools.org/latest/gmt.html#w-full.

Returns

track (pandas.DataFrame or None) – Return type depends on whether the outfile parameter is set:

• pandas.DataFrame table with (x, y, …, newcolname) if outfile is not set

• None if outfile is set (track output will be stored in file set by outfile)

Example

>>> import pygmt
>>> # Load a grid of @earth_relief_30m data, with an x-range of -118 to
>>> # -107, and a y-range of -49 to -42
>>> grid = pygmt.datasets.load_earth_relief(
...     resolution="30m", region=[-118, -107, -49, -42]
... )
>>> # Load a pandas dataframe with ocean ridge points
>>> points = pygmt.datasets.load_sample_data(name="ocean_ridge_points")
>>> # Create a pandas dataframe from an input grid and set of points
>>> # The output dataframe adds a column named "bathymetry"
>>> output_dataframe = pygmt.grdtrack(
...     points=points, grid=grid, newcolname="bathymetry"
... )

Examples using pygmt.grdtrack

Sampling along tracks