Provided by: gmt-common_5.4.5+dfsg-2_all bug

NAME

       gmtspatial - Do geospatial operations on lines and polygons

SYNOPSIS

       gmtspatial     [     table     ]     [      -A[amin_dist][unit]]     [      -C     ]     [
       -D[+ffile][+aamax][+ddmax][+c|Ccmax][+sfact] ] [  -E+|- ]  [   -F[l]  ]  [   -I[e|i]  ]  [
       -Npfile[+a][+pstart][+r][+z]  ]  [   -Q[[-|+]unit][+cmin[/max]][+h][+l][+p][+s[a|d]]  ]  [
       -Rregion ] [  -Si|u|s|j ] [  -T[clippolygon] ] [  -V[level] ] [ -bbinary ] [ -dnodata ]  [
       -eregexp ] [ -fflags ] [ -ggaps ] [ -hheaders ] [ -iflags ] [ -oflags ] [ -:[i|o] ]

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

DESCRIPTION

       gmtspatial  reads  one  or more data files (which may be multisegment files) that contains
       closed polygons and operates of these polygons in the specified  way.  Operations  include
       area calculation, handedness reversals, and polygon intersections.

REQUIRED ARGUMENTS

       None.

OPTIONAL ARGUMENTS

       table  One  or  more  ASCII (or binary, see -bi[ncols][type]) data table file(s) holding a
              number of data columns. If no tables are given then we read from standard input.

       -A[amin_dist][unit]
              Perform spatial nearest neighbor (NN) analysis: Determine the nearest  neighbor  of
              each point and report the NN distances and the point IDs involved in each pair (IDs
              are the input record numbers starting at 0).  Use -Aa to decimate  a  data  set  so
              that  no  NN  distance is lower than the threshold min_dist.  In this case we write
              out the (possibly averaged) coordinates and the updated NN distances and point IDs.
              A negative point number means the original point was replaced by a weighted average
              (the absolute ID value gives the ID of the first original point ID to  be  included
              in  the  average.).   Note: The input data are assumed to contain (lon, lat) or (x,
              y), optionally followed by a z and a weight [1]  column.   We  compute  a  weighted
              average of the location and z (if present).

       -C     Clips  polygons to the map region, including map boundary to the polygon as needed.
              The result is a closed polygon (see -T for truncation instead). Requires -R.

       -D[+ffile][+aamax][+ddmax][+c|Ccmax][+sfact]
              Check for duplicates among the input lines or polygons, or, if file  is  given  via
              +f,  check  if  the  input  features  already  exist among the features in file. We
              consider the cases of exact (same number and coordinates) and  approximate  matches
              (average distance between nearest points of two features is less than a threshold).
              We also consider that some features may have been reversed. Features are considered
              approximate matches if their minimum distance is less than dmax [0] (see UNITS) and
              their closeness (defined as the ratio between  the  average  distance  between  the
              features  divided  by  their  average  length)  is  less than cmax [0.01]. For each
              duplicate found, the output record begins with the single letter Y (exact match) or
              ~ (approximate match). If the two matching segments differ in length by more than a
              factor of 2 then we consider the duplicate to be either a subset (-) or a  superset
              (+).  Finally,  we  also note if two lines are the result of splitting a continuous
              line across the Dateline  (|).   For  polygons  we  also  consider  the  fractional
              difference  in  areas; duplicates must differ by less than amax [0.01]. By default,
              we compute the mean line separation. Use +Ccmin to instead compute the median  line
              separation  and therefore a robust closeness value. Also by default we consider all
              distances between points on one line and another. Append +p to limit the comparison
              to  points  that  project  perpendicularly to points on the other line (and not its
              extension).

       -E+|- ]
              Reset the handedness of all polygons to match the given + (counter-clockwise) or  -
              (clockwise). Implies -Q+.

       -F[l]  Force  input  data to become polygons on output, i.e., close them explicitly if not
              already closed.  Optionally, append l to force line geometry.

       -I[e|i]
              Determine the intersection locations between all pairs of polygons.   Append  i  to
              only  compute  internal  (i.e., self-intersecting polygons) crossovers or e to only
              compute external (i.e., between paris of polygons) crossovers [Default is both].

       -Npfile[+a][+pstart][+r][+z]
              Determine if one (or all, with +a) points of each feature in  the  input  data  are
              inside any of the polygons given in the pfile. If inside, then report which polygon
              it is; the polygon ID is either taken from the aspatial value assigned  to  Z,  the
              segment  header  (first  -Z,  then  -L  are scanned), or it is assigned the running
              number that is initialized to start [0]. By default  the  input  segment  that  are
              found  to  be inside a polygon are written to stdout with the polygon ID encoded in
              the segment header as -ZID. Alternatively, append +r to just report  which  polygon
              contains  a  feature or +z to have the IDs added as an extra data column on output.
              Segments that fail to be inside a polygon are not written out.  If  more  than  one
              polygon contains the same segment we skip the second (and further) scenario.

       -Q[[-|+]unit][+cmin[/max]][+h][+l][+p][+s[a|d]]
              Measure the area of all polygons or length of line segments. Use -Q+h to append the
              area to each polygons segment header [Default simply writes the  area  to  stdout].
              For  polygons  we also compute the centroid location while for line data we compute
              the mid-point (half-length) position. Append a distance unit  to  select  the  unit
              used  (see  UNITS).  Note  that  the  area  will  depend  on the current setting of
              PROJ_ELLIPSOID; this should be a recent ellipsoid  to  get  accurate  results.  The
              centroid  is  computed  using  the  mean of the 3-D Cartesian vectors making up the
              polygon vertices, while the area is obtained via  an  equal-area  projection.   For
              line  lengths  you  may  prepend  -|+ to the unit and the calculation will use Flat
              Earth or Geodesic algorithms, respectively [Default  is  great  circle  distances].
              Normally,  all  input  segments  will be be reflected on output.  Use c to restrict
              processing to those whose length (or area for polygons) fall inside  the  specified
              range  set by min and max.  If max is not set it defaults to infinity.  To sort the
              segments based on their lengths or area, use s and append a for ascending and d for
              descending  order  [ascending].   By  default,  we consider open polygons as lines.
              Append +p to close open polygons and thus consider all input as polygons, or append
              +l to consider all input as lines, even if closed.

       -Rwest/east/south/north[/zmin/zmax][+r][+uunit]
              west,  east,  south,  and north specify the region of interest, and you may specify
              them in decimal degrees or in [±]dd:mm[:ss.xxx][W|E|S|N] format Append +r if  lower
              left  and  upper  right  map  coordinates  are  given  instead  of w/e/s/n. The two
              shorthands -Rg and -Rd stand for global domain (0/360 and  -180/+180  in  longitude
              respectively,  with  -90/+90  in  latitude).  Alternatively for grid creation, give
              Rcodelon/lat/nx/ny, where code is a 2-character combination of L, C, R  (for  left,
              center,  or right) and T, M, B for top, middle, or bottom. e.g., BL for lower left.
              This indicates which point on a rectangular region the  lon/lat  coordinate  refers
              to,  and  the grid dimensions nx and ny with grid spacings via -I is used to create
              the corresponding region.  Alternatively, specify the name of an existing grid file
              and  the  -R  settings  (and grid spacing, if applicable) are copied from the grid.
              Appending +uunit expects projected (Cartesian) coordinates compatible  with  chosen
              -J and we inversely project to determine actual rectangular geographic region.  For
              perspective view (-p), optionally append /zmin/zmax.  In case of  perspective  view
              (-p),  a z-range (zmin, zmax) can be appended to indicate the third dimension. This
              needs to be done only when using the -Jz option, not when using only the -p option.
              In  the  latter  case  a  perspective  view  of the plane is plotted, with no third
              dimension. Clips polygons to the map region, including map boundary to the  polygon
              as needed. The result is a closed polygon.

       -Si|j|s|u
              Spatial  processing  of polygons. Choose from -Si which returns the intersection of
              polygons (closed), -Su which returns the union of polygons (closed), -Ss which will
              split  polygons  that  straddle the Dateline, and -Sj which will join polygons that
              were split by the Dateline.  Note: Only -Ss has been implemented.

       -T[clippolygon]
              Truncate polygons against the specified polygon given, possibly resulting  in  open
              polygons.  If no argument is given to -T we create a clipping polygon from -R which
              then is required. Note that when the -R clipping is in effect we will also look for
              polygons of length 4 or 5 that exactly match the -R clipping polygon.

       -V[level] (more ...)
              Select verbosity level [c].

       -bi[ncols][t] (more ...)
              Select native binary input. [Default is 2 input columns].

       -bo[ncols][type] (more ...)
              Select native binary output. [Default is same as input].

       -d[i|o]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.

       -f[i|o]colinfo (more ...)
              Specify data types of input and/or output columns.

       -g[a]x|y|d|X|Y|D|[col]z[+|-]gap[u] (more ...)
              Determine data gaps and line breaks.

       -h[i|o][n][+c][+d][+rremark][+rtitle] (more ...)
              Skip or produce header record(s).

       -icols[+l][+sscale][+ooffset][,...] (more ...)
              Select input columns and transformations (0 is first column).

       -ocols[,...] (more ...)
              Select output columns (0 is first column).

       -:[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 exits (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 exits.

       -? or no arguments
              Print  a  complete  usage (help) message, including the explanation of all options,
              then exits.

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], 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. Prepend - to a distance (or the unit is no distance is
       given) to perform "Flat Earth" calculations (quicker but less accurate) or  prepend  +  to
       perform exact geodesic calculations (slower but more accurate).

ASCII FORMAT PRECISION

       The  ASCII  output formats of numerical data are controlled by parameters in your gmt.conf
       file. Longitude and latitude are formatted according to FORMAT_GEO_OUT, absolute  time  is
       under  the control of FORMAT_DATE_OUT and FORMAT_CLOCK_OUT, whereas general floating point
       values are formatted according to FORMAT_FLOAT_OUT. Be aware that the format in effect can
       lead  to loss of precision in ASCII output, which can lead to various problems downstream.
       If you find the output is not written with enough precision, consider switching to  binary
       output (-bo if available) or specify more decimals using the FORMAT_FLOAT_OUT setting.

EXAMPLE

       To turn all lines in the multisegment file lines.txt into closed polygons, run

              gmt spatial lines.txt -F > polygons.txt

       To compute the area of all geographic polygons in the multisegment file polygons.txt, run

              gmt spatial polygons.txt -Q > areas.txt

       Same  data,  but  now orient all polygons to go counter-clockwise and write their areas to
       the segment headers, run

              gmt spatial polygons.txt -Q+h -E+ > areas.txt

       To determine the areas of all the polygon segments in the file janmayen_land_full.txt, add
       this  information  to  the  segment headers, sort the segments from largest to smallest in
       area but only keep polygons with area larger than 1000 sq. meters, run

              gmt spatial -Qe+h+p+c1000+sd -V janmayen_land_full.txt > largest_pols.txt

       To determine the intersections between the polygons A.txt and B.txt, run

              gmt spatial A.txt B.txt -Ie > crossovers.txt

       To truncate polygons A.txt against polygon B.txt, resulting in an open line segment, run

              gmt gmtspatial A.txt -TB.txt > line.txt

NOTES

       OGR/GMT files are considered complete datasets and thus you cannot specify more  than  one
       at  a  given  time.  This  causes problems if you want to examine the intersections of two
       OGR/GMT files.  The solution is to convert them to regular  datasets  via  gmtconvert  and
       then run gmtspatial on the converted files.

SEE ALSO

       gmt, gmtconvert, gmtselect, gmtsimplify

COPYRIGHT

       2019, P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis, and F. Wobbe