some issues with -F option on gmtselect
|Target version:||Candidate for next minor release|
The following have been tested using gmtselect(core) 5.4.0_r17525M [64-bit]
There seems to be some issues with the -F option on gmtselect.
Attached are two polygon files ("lat lon" in each record) both describing a circle on a sphere of 15 degree radius with point separation of 0.5 degree. On of the polygons is centred on lat/lon 75/0 (hence have a point at the north pole) the other is centred on lat/lon 76/0 (hence encloses the north pole)
Consider the polygon centerd at 76/0 as the polygon to select points inside/outside and the points of the polygon centered at 75/0 as the points to select from.
To get the points inside or on the polygon we run the command:
gmtselect -: 75_0_15_05.poly -F76_0_15_05.poly > tmp_inpoints.latlon
which results in 250 out of the 720 points being selected
To get hte points outside the polygon we instead run the command:
gmtselect -: 75_0_15_05.poly -If -F76_0_15_05.poly > tmp_outpoints.latlon
which results in 420 out of the 720 points being selected
To visualize the result we run the commands:
gmtset MAP_GRID_CROSS_SIZE_PRIMARY 0.1i MAP_FRAME_TYPE FANCY PROJ_LENGTH_UNIT cm PS_MEDIA a4
pscoast -Rg -JG0/75/320000/0/0/0/0/0/15 -Dc -A1000 -Glightbrown -Wthinnest -P -Slightblue -K > fig.ps
psxy -: 76_0_15_05.poly -R -J -W1,black -O -K >> fig.ps
psxy -: 76_0_15_05.poly -R -J -Sx0.05 -W4,black -O -K >> fig.ps
psxy -: tmp_inpoints.latlon -R -J -Sc0.05 -W1,blue -O -K >> fig.ps
psxy -: tmp_outpoints.latlon -R -J -Sc0.05 -W1,red -O -K >> fig.ps
psbasemap -R -J -Bagf -BWSen -O >> fig.ps
The result can be seen in the attached figure where points inside are marked in read and outside in blue, clearly there are points the are erroneously considered to be inside the polygon when actually outside (see red points close to or at latitude 60).
Note, if we swap which file is considered to be the polygon and which is considered to be the set of points all works as should. I suspect the problem resides in the fact that in the above case the polygon encloses the north pole.
Further if the same points are used for the polygon as for the points to select from gmtselect hangs indefinitely (well at least long enough for me to get tired and kill it) i.e. running the command:
gmtselect -: 75_0_15_05.poly -F75_0_15_05.poly > tmp.latlon
Now, while this may seem a stupid thing to do strictly gmtselect should not hang but return all the points as these are exactly on the polygon
Ok, sorry missed that. Think I assumed that gmtselect always assumed geographical coordinates as it is stated under "Description" that "gmtselect is a filter that reads (longitude, latitude) positions from the first 2 columns of infiles" then again under the -F option a spherical polygon is mentioned to require specific considerations (with regards to separation between points) so perhaps this should have been understood to also imply the need for -fg.
However, to aid future users perhaps one could mention under the description of the -F option that -fg may be required, or better yet have gmt always assume geographical coordinates as this seems natural and opt for -fp when using cartesian coordinates (personally I prefer this).
But, even with the -fg option gmtselect seems to hang if the points to select from are the same as the points in the polygon, i.e. using the command:
gmtselect -: 75_0_15_05.poly -F75_0_15_05.poly
Further, if a point of the polygon is located at the north pole (probably also the south pole, but haven't tested) gmtselect fails to identify a point to be on the polygon unless the longitude is identical to that in the input file, e.g.:
$~> echo "90 30" | gmtselect -: -F75_0_15_05.poly -fg
$~> echo "90 90" | gmtselect -: -F75_0_15_05.poly -fg
$~> echo "90 180" | gmtselect -: -F75_0_15_05.poly -fg
$~> grep 90.000 75_0_15_05.poly
$~> echo "90 0" | gmtselect -: -F75_0_15_05.poly -fg
I guess these are minor issues that may not arise very frequently