1GMTSPATIAL(1) GMT GMTSPATIAL(1)
2
6 gmtspatial - Do geospatial operations on lines and polygons
7
9 gmtspatial [ table ] [ -A[amin_dist][unit]] [ -C ] [
10 -D[+ffile][+aamax][+ddmax][+c|Ccmax][+sfact] ] [ -E+|- ] [ -F[l] ] [
11 -I[e|i] ] [ -Npfile[+a][+pstart][+r][+z] ] [
12 -Q[[-|+]unit][+cmin[/max]][+h][+l][+p][+s[a|d]] ] [ -Rregion ] [
13 -Si|u|s|j ] [ -T[clippolygon] ] [ -V[level] ] [ -bbinary ] [ -dnodata
14 ] [ -eregexp ] [ -fflags ] [ -ggaps ] [ -hheaders ] [ -iflags ] [
15 -oflags ] [ -:[i|o] ]
16 Note: No space is allowed between the option flag and the associated
18 arguments.
19
21 gmtspatial reads one or more data files (which may be multisegment
22 files) that contains closed polygons and operates of these polygons in
23 the specified way. Operations include area calculation, handedness
24 reversals, and polygon intersections.
25
27 None.
28
30 table One or more ASCII (or binary, see -bi[ncols][type]) data table
31 file(s) holding a number of data columns. If no tables are given
32 then we read from standard input.
33 -A[amin_dist][unit]
35 Perform spatial nearest neighbor (NN) analysis: Determine the
36 nearest neighbor of each point and report the NN distances and
37 the point IDs involved in each pair (IDs are the input record
38 numbers starting at 0). Use -Aa to decimate a data set so that
39 no NN distance is lower than the threshold min_dist. In this
40 case we write out the (possibly averaged) coordinates and the
41 updated NN distances and point IDs. A negative point number
42 means the original point was replaced by a weighted average (the
43 absolute ID value gives the ID of the first original point ID to
44 be included in the average.). Note: The input data are assumed
45 to contain (lon, lat) or (x, y), optionally followed by a z and
46 a weight [1] column. We compute a weighted average of the loca‐
47 tion and z (if present).
48 -C Clips polygons to the map region, including map boundary to the
50 polygon as needed. The result is a closed polygon (see -T for
51 truncation instead). Requires -R.
52 -D[+ffile][+aamax][+ddmax][+c|Ccmax][+sfact]
54 Check for duplicates among the input lines or polygons, or, if
55 file is given via +f, check if the input features already exist
56 among the features in file. We consider the cases of exact (same
57 number and coordinates) and approximate matches (average dis‐
58 tance between nearest points of two features is less than a
59 threshold). We also consider that some features may have been
60 reversed. Features are considered approximate matches if their
61 minimum distance is less than dmax [0] (see UNITS) and their
62 closeness (defined as the ratio between the average distance
63 between the features divided by their average length) is less
64 than cmax [0.01]. For each duplicate found, the output record
65 begins with the single letter Y (exact match) or ~ (approximate
66 match). If the two matching segments differ in length by more
67 than a factor of 2 then we consider the duplicate to be either a
68 subset (-) or a superset (+). Finally, we also note if two lines
69 are the result of splitting a continuous line across the Date‐
70 line (|). For polygons we also consider the fractional differ‐
71 ence in areas; duplicates must differ by less than amax [0.01].
72 By default, we compute the mean line separation. Use +Ccmin to
73 instead compute the median line separation and therefore a
74 robust closeness value. Also by default we consider all dis‐
75 tances between points on one line and another. Append +p to
76 limit the comparison to points that project perpendicularly to
77 points on the other line (and not its extension).
78 -E+|- ]
80 Reset the handedness of all polygons to match the given +
81 (counter-clockwise) or - (clockwise). Implies -Q+.
82 -F[l] Force input data to become polygons on output, i.e., close them
84 explicitly if not already closed. Optionally, append l to force
85 line geometry.
86 -I[e|i]
88 Determine the intersection locations between all pairs of poly‐
89 gons. Append i to only compute internal (i.e., self-intersect‐
90 ing polygons) crossovers or e to only compute external (i.e.,
91 between paris of polygons) crossovers [Default is both].
92 -Npfile[+a][+pstart][+r][+z]
94 Determine if one (or all, with +a) points of each feature in the
95 input data are inside any of the polygons given in the pfile. If
96 inside, then report which polygon it is; the polygon ID is
97 either taken from the aspatial value assigned to Z, the segment
98 header (first -Z, then -L are scanned), or it is assigned the
99 running number that is initialized to start [0]. By default the
100 input segment that are found to be inside a polygon are written
101 to stdout with the polygon ID encoded in the segment header as
102 -ZID. Alternatively, append +r to just report which polygon con‐
103 tains a feature or +z to have the IDs added as an extra data
104 column on output. Segments that fail to be inside a polygon are
105 not written out. If more than one polygon contains the same seg‐
106 ment we skip the second (and further) scenario.
107 -Q[[-|+]unit][+cmin[/max]][+h][+l][+p][+s[a|d]]
109 Measure the area of all polygons or length of line segments. Use
110 -Q+h to append the area to each polygons segment header [Default
111 simply writes the area to stdout]. For polygons we also compute
112 the centroid location while for line data we compute the
113 mid-point (half-length) position. Append a distance unit to
114 select the unit used (see UNITS). Note that the area will depend
115 on the current setting of PROJ_ELLIPSOID; this should be a
116 recent ellipsoid to get accurate results. The centroid is com‐
117 puted using the mean of the 3-D Cartesian vectors making up the
118 polygon vertices, while the area is obtained via an equal-area
119 projection. For line lengths you may prepend -|+ to the unit
120 and the calculation will use Flat Earth or Geodesic algorithms,
121 respectively [Default is great circle distances]. Normally, all
122 input segments will be be reflected on output. Use c to
123 restrict processing to those whose length (or area for polygons)
124 fall inside the specified range set by min and max. If max is
125 not set it defaults to infinity. To sort the segments based on
126 their lengths or area, use s and append a for ascending and d
127 for descending order [ascending]. By default, we consider open
128 polygons as lines. Append +p to close open polygons and thus
129 consider all input as polygons, or append +l to consider all
130 input as lines, even if closed.
131 -Rwest/east/south/north[/zmin/zmax][+r][+uunit]
133 west, east, south, and north specify the region of interest, and
134 you may specify them in decimal degrees or in
135 [±]dd:mm[:ss.xxx][W|E|S|N] format Append +r if lower left and
136 upper right map coordinates are given instead of w/e/s/n. The
137 two shorthands -Rg and -Rd stand for global domain (0/360 and
138 -180/+180 in longitude respectively, with -90/+90 in latitude).
139 Alternatively for grid creation, give Rcodelon/lat/nx/ny, where
140 code is a 2-character combination of L, C, R (for left, center,
141 or right) and T, M, B for top, middle, or bottom. e.g., BL for
142 lower left. This indicates which point on a rectangular region
143 the lon/lat coordinate refers to, and the grid dimensions nx and
144 ny with grid spacings via -I is used to create the corresponding
145 region. Alternatively, specify the name of an existing grid
146 file and the -R settings (and grid spacing, if applicable) are
147 copied from the grid. Appending +uunit expects projected (Carte‐
148 sian) coordinates compatible with chosen -J and we inversely
149 project to determine actual rectangular geographic region. For
150 perspective view (-p), optionally append /zmin/zmax. In case of
151 perspective view (-p), a z-range (zmin, zmax) can be appended to
152 indicate the third dimension. This needs to be done only when
153 using the -Jz option, not when using only the -p option. In the
154 latter case a perspective view of the plane is plotted, with no
155 third dimension. Clips polygons to the map region, including map
156 boundary to the polygon as needed. The result is a closed poly‐
157 gon.
158 -Si|j|s|u
160 Spatial processing of polygons. Choose from -Si which returns
161 the intersection of polygons (closed), -Su which returns the
162 union of polygons (closed), -Ss which will split polygons that
163 straddle the Dateline, and -Sj which will join polygons that
164 were split by the Dateline. Note: Only -Ss has been imple‐
165 mented.
166 -T[clippolygon]
168 Truncate polygons against the specified polygon given, possibly
169 resulting in open polygons. If no argument is given to -T we
170 create a clipping polygon from -R which then is required. Note
171 that when the -R clipping is in effect we will also look for
172 polygons of length 4 or 5 that exactly match the -R clipping
173 polygon.
174 -V[level] (more ...)
176 Select verbosity level [c].
177 -bi[ncols][t] (more ...)
179 Select native binary input. [Default is 2 input columns].
180 -bo[ncols][type] (more ...)
182 Select native binary output. [Default is same as input].
183 -d[i|o]nodata (more ...)
185 Replace input columns that equal nodata with NaN and do the
186 reverse on output.
187 -e[~]"pattern" | -e[~]/regexp/[i] (more ...)
189 Only accept data records that match the given pattern.
190 -f[i|o]colinfo (more ...)
192 Specify data types of input and/or output columns.
193 -g[a]x|y|d|X|Y|D|[col]z[+|-]gap[u] (more ...)
195 Determine data gaps and line breaks.
196 -h[i|o][n][+c][+d][+rremark][+rtitle] (more ...)
198 Skip or produce header record(s).
199 -icols[+l][+sscale][+ooffset][,...] (more ...)
201 Select input columns and transformations (0 is first column).
202 -ocols[,...] (more ...)
204 Select output columns (0 is first column).
205 -:[i|o] (more ...)
207 Swap 1st and 2nd column on input and/or output.
208 -^ or just -
210 Print a short message about the syntax of the command, then
211 exits (NOTE: on Windows just use -).
212 -+ or just +
214 Print an extensive usage (help) message, including the explana‐
215 tion of any module-specific option (but not the GMT common
216 options), then exits.
217 -? or no arguments
219 Print a complete usage (help) message, including the explanation
220 of all options, then exits.
221
223 For map distance unit, append unit d for arc degree, m for arc minute,
224 and s for arc second, or e for meter [Default], f for foot, k for km, M
225 for statute mile, n for nautical mile, and u for US survey foot. By
226 default we compute such distances using a spherical approximation with
227 great circles. Prepend - to a distance (or the unit is no distance is
228 given) to perform "Flat Earth" calculations (quicker but less accurate)
229 or prepend + to perform exact geodesic calculations (slower but more
230 accurate).
231
233 The ASCII output formats of numerical data are controlled by parameters
234 in your gmt.conf file. Longitude and latitude are formatted according
235 to FORMAT_GEO_OUT, absolute time is under the control of FOR‐
236 MAT_DATE_OUT and FORMAT_CLOCK_OUT, whereas general floating point val‐
237 ues are formatted according to FORMAT_FLOAT_OUT. Be aware that the for‐
238 mat in effect can lead to loss of precision in ASCII output, which can
239 lead to various problems downstream. If you find the output is not
240 written with enough precision, consider switching to binary output (-bo
241 if available) or specify more decimals using the FORMAT_FLOAT_OUT set‐
242 ting.
243
245 To turn all lines in the multisegment file lines.txt into closed poly‐
246 gons, run
247 gmt spatial lines.txt -F > polygons.txt
249 To compute the area of all geographic polygons in the multisegment file
251 polygons.txt, run
252 gmt spatial polygons.txt -Q > areas.txt
254 Same data, but now orient all polygons to go counter-clockwise and
256 write their areas to the segment headers, run
257 gmt spatial polygons.txt -Q+h -E+ > areas.txt
259 To determine the areas of all the polygon segments in the file jan‐
261 mayen_land_full.txt, add this information to the segment headers, sort
262 the segments from largest to smallest in area but only keep polygons
263 with area larger than 1000 sq. meters, run
264 gmt spatial -Qe+h+p+c1000+sd -V janmayen_land_full.txt > largest_pols.txt
266 To determine the intersections between the polygons A.txt and B.txt,
268 run
269 gmt spatial A.txt B.txt -Ie > crossovers.txt
271 To truncate polygons A.txt against polygon B.txt, resulting in an open
273 line segment, run
274 gmt gmtspatial A.txt -TB.txt > line.txt
276
278 OGR/GMT files are considered complete datasets and thus you cannot
279 specify more than one at a given time. This causes problems if you want
280 to examine the intersections of two OGR/GMT files. The solution is to
281 convert them to regular datasets via gmtconvert and then run gmtspatial
282 on the converted files.
283
285 gmt, gmtconvert, gmtselect, gmtsimplify
286
288 2019, P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis, and F. Wobbe
2895.4.5 Feb 24, 2019 GMTSPATIAL(1)