Augment gmtspatial with -N+i

[PaulWessel]

The -N+i setting will treat each input point individually and determine which polygon it falls into (if any). If a polygon ID is found then we write that point's full record out and add a final column with the ID.

Test example with some fake points and polygons::

cat << EOF > cloud.txt
138 15
139 15
140 15
141 15
142 15
143 15
138 12
139 12
141 18
142 18
150 20
EOF
cat << EOF > pols.txt
> -Z999
140 10
140 20
150 20
150 10
140 10
> -Z888
130 10
130 20
140 20
140 10
130 10
EOF

Then running the command:

gmt spatial test_points.txt -Npols.txt+i
141	15	999
142	15	999
143	15	999
141	18	999
142	18	999
138	15	888
139	15	888
138	15	888
138	12	888
139	12	888

So 5 points each in polygon IDs 888 and 999, plus one point outside all and hence skipped.
Note: I have labelled it WIP and would appreciate some more testing before we approve/merge.

[joa-quim]

Gave the same result as my (very) simple test for the GMT.jl inwhichpolygon() function.

[joa-quim]

Wait, BUG. If headers have no -Z or -L

gmt spatial cloud.txt -Npols.txt+i
gmtspatial [ERROR]: No polygon ID found; ID set to NaN
141     15      -1
142     15      -1
143     15      -1
141     18      -1
142     18      -1
gmtspatial [ERROR]: No polygon ID found; ID set to NaN
138     15      -1
139     15      -1
138     12      -1
139     12      -1
138     15      -1

And I would really love to have a way to get ID's start counting at 1, not 0.

[PaulWessel]

Well, see +p which we then need to allow

[PaulWessel]

On bus so hard to check

[joa-quim]

Yes, +p can be the solution.