grdmask

grdmask(cmd0::String="", arg1=nothing, kwargs...)

keywords: GMT, Julia, masks, grid

Create mask grid from polygons or point coverage.

Description

grdmask can operate in two different modes. 1. It reads one or more pathfiles that each define a closed polygon. The nodes defined by the specified region and lattice spacing will be set equal to one of three possible values depending on whether the node is outside, on the polygon perimeter, or inside the polygon, with the assigned z value selected via out_edge_in. If multiple polygons overlap the same nodes then the polygon selected depends on the clobber selection. The resulting mask may be used in subsequent operations involving grdmath to mask out data from polygonal areas. 2. The pathfiles simply represent data point locations and the mask is set to the inside or outside value depending on whether a node is within a maximum distance from the nearest data point. If the distance specified is zero then only the nodes nearest each data point are considered "inside".

Required Arguments

• table
  One or more data tables holding a number of data columns.

• I or inc or increment or spacing : – inc=x_inc | inc=(xinc, yinc) | inc="xinc[+e|n][/yinc[+e|n]]"
  Specify the grid increments or the block sizes. More at spacing

• R or region or limits : – limits=(xmin, xmax, ymin, ymax) | limits=(BB=(xmin, xmax, ymin, ymax),) | limits=(LLUR=(xmin, xmax, ymin, ymax),units="unit") | ...more
  Specify the region of interest. More at limits. For perspective view view, optionally add zmin,zmax. This option may be used to indicate the range used for the 3-D axes. You may ask for a larger w/e/s/n region to have more room between the image and the axes.

Optional Arguments

• A or steps : – steps=true | steps=:meridian|:parallel|:x|:y|:r|:theta
  By default, geographic line segments (as indicated for example by the colinfo option) are drawn as great circle arcs. To draw them as straight lines, use the steps=true. Alternatively, use steps=:meridian to draw the line by first following a meridian, then a parallel. Or append steps=:parallel to start following a parallel, then a meridian. (This can be practical to draw a line along parallels, for example). For Cartesian data, points are simply connected, unless you use steps=:x or steps=:y to draw stair-case curves that whose first move is along x or y, respectively. If your Cartesian data are polar (theta, r), use steps=:t or steps=:r to construct stair-case paths whose first move is along theta or r, respectively.

• C or clobber : – clobber=:f|l|o|u
  Clobber mode: Selects the polygon whose z-value will determine the grid nodes. Choose from the following modes: clobber=:f for the first polygon to overlap a node; o for the last polygon to overlap a node; l for the polygon with the lowest z-value, and u for the polygon with the uppermost z-value [Default is o]. Note: Does not apply to search_radius. For polygon z-values, see outedgein.

• G or save or outgrid or outfile : – outgrid=[=ID][+ddivisor][+ninvalid][+ooffset|a][+sscale|a][:driver[dataType][+coptions]]
  Give the name of the output grid file. Optionally, append =ID for writing a specific file format (See full description). 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].

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

Note2: This is optional and to be used only when saving the result directly on disk. Otherwise, just use the G = modulename(...) form.

• N or out_edge_in : – out_edge_in=values | out_edge_in="z|Z|p|Pvalues"
  Sets the out/edge/in that will be assigned to nodes that are outside the polygons, on the edge, or inside. Values can be any number, including the textstring NaN [Default is 0/0/1]. Optionally, use Nz to set polygon insides to the z-value obtained from the data (either segment header Zzval, Lheader or via aZ=name); use NZ to consider the polygon boundary as part of the inside. Alternatively, use Np to use a running number as polygon ID; optionally append start of the sequence [0]. Here, NP includes the polygon perimeter as inside. Note: out_edge_in="z|Z|p|P" cannot be used in conjunction with seach_radius; they also all optionally accept /out [0].

• S or search_radius : – search_radius=rad | search_radius=(xlim,ylim)
  Set nodes to inside, on edge, or outside depending on their distance to the nearest data point. Nodes within radius [0] from the nearest data point are considered inside; append a distance unit (see Units). If radius is given as z then we instead read individual radii from the 3rd input column. Unless Cartesian data, specify the unit of these radii by appending it after -Sz. If search_radius is not set then we consider the input data to define one or more closed polygon(s) instead. For Cartesian data with different units you can instead use search_radius=(xlim,ylim) which will perform a rectangular search where all nodes within ±\xlim and ±\ylim of a data point will be considered inside. One can also achieve the rectangular selection effect by using the search_radius=ncellsc form. Here *ncells* means the number of cells around each data point. As an example, search_radius="0c" means that only the cell where point lies is masked, search_radius="1c" masks one cell beyond that (i.e. makes a 3x3 neighborhood), and so on.

• V or verbose : – verbose=true | verbose=level
  Select verbosity level. More at verbose

• a or aspatial : – aspatial=??
  Control how aspatial data are handled in GMT during input and output. More at

• bi or binary_in : – binary_in=??
  Select native binary format for primary table input. More at

• di or nodata_in : – nodata_in=??
  Substitute specific values with NaN. More at

• e or pattern : – pattern=??
  Only accept ASCII data records that contain the specified pattern. More at

• f or colinfo : – colinfo=??
  Specify the data types of input and/or output columns (time or geographical data). More at

• g or gap : – gap=??
  Examine the spacing between consecutive data points in order to impose breaks in the line. More at

• h or header : – header=??
  Specify that input and/or output file(s) have n header records. More at

• i or incol or incols : – incol=col_num | incol="opts"
  Select input columns and transformations (0 is first column, t is trailing text, append word to read one word only). More at incol

• j or metric or spherical_dist or spherical : – metric=greatcirc or spherical=:flat or spherical=:ellipsoidal
  Determine how spherical distances are calculated in modules that support this [Default is spherical=:greatcirc]. GMT has different ways to compute distances on planetary bodies:

- spherical=:greatcirc to perform great circle distance calculations, with parameters such as distance increments or radii compared against calculated great circle distances [Default is spherical=:greatcirc]. - spherical=:flat to select Flat Earth mode, which gives a more approximate but faster result. - spherical=:ellipsoidal to select ellipsoidal (or geodesic) mode for the highest precision and slowest calculation time. Note: 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).

• n or interp or interpol : – interp=params
  Select interpolation mode for grids. More at interp

• q or inrows : – inrows=??
  Select specific data rows to be read and/or written. More at

• r or reg or registration : – reg=:p | reg=:g
  Select gridline or pixel node registration. Used only when output is a grid. More at

• w or wrap or cyclic : – wrap=??
  Convert input records to a cyclical coordinate. More at

• x or cores : – cores=N
  Limit the number of cores used when making the individual frames. By default we try to use all available cores. Append N to only use n cores (if too large it will be truncated to the maximum cores available). Finally, give a negative N to select (all - N) cores (or at least 1 if N equals or exceeds all). The parallel processing does not depend on OpenMP.

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

Inside/outside Status

To determine if a point is inside, outside, or exactly on the boundary of a polygon we need to balance the complexity (and execution time) of the algorithm with the type of data and shape of the polygons. For any Cartesian data we use a non-zero winding algorithm, which is quite fast. For geographic data we will also use this algorithm as long as (1) the polygons do not include a geographic pole, and (2) the longitude extent of the polygons is less than 360. If this is the situation we also carefully adjust the test point longitude for any 360 degree offsets, if appropriate. Otherwise, we employ a full spherical ray-shooting method to determine a points status.

Notes

A grid produced by grdmask is a categorical dataset. As such, one has to be careful not to interpolate it with standard methods, such as splines. However, if you make a map of this grid using a map projection the grid will be reprojected to yield a rectangular matrix in the projected coordinates. This interpolation is done using splines by default and thus may yield artifacts in your map. We recommend you use grdimageinterp=:n to instead use a nearest neighbor interpolation for such cases.

Save storage space

Since most uses of grdmask revolves around creating mask grids that hold just a few integer values (and perhaps NaN), we choose to write them to disk as byte grids by appending the suffix =nb to the desired grid filename. Some situations may store integers that exceed the range available in a byte and for those we specify a short integer grid with =ns. For larger integers you may consider =ni, otherwise use the default float grid format.

Examples

To set all nodes inside and on the polygons coastline_*.xy to 0, and outside points to 1, do

G = grdmask("coastline_*.xy", region=(-60,-40,-40,-30), inc="5m", out_edge_in=(1,0,0))

To set nodes within 50 km of data points to 1 and other nodes to NaN, do

G = grdmask("data.xyz", region=(-60,-40,-40,-30), inc="5m", out_edge_in=(NaN,1,1), search_radius="50k")

To assign polygon IDs to the gridnodes using the insides of the polygons in plates.gmt, based on the attribute POL_ID, do

G = grdmask("plates.gmt", region=(-40,-40,-40,-40), inc="2m", out_edge_in=:z, aZ=POL_ID)

Same exercise, but instead compute running polygon IDs starting at 100, do

G = grdmask("plates.gmt", region=(-40,-40,-40,-40), inc="2m", out_edge_in="p100")

See Also

grdlandmask, grdclip, mask, clip