Bug #1085

some issues with -F option on gmtselect

Added by Peter 7 months ago. Updated 6 months ago.

Status:ClosedStart date:2017-04-28
Priority:NormalDue date:
Assignee:Paul% Done:


Target version:Candidate for next minor release
Affected version:5.x-svn Platform:



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

75_0_15_05.poly - polygon, circle centered at lat/lon 75/0, radius 15 degree, point separation 0.5 degree (14.1 KB) Peter, 2017-04-28 10:54

76_0_15_05.poly - polygon, circle centered at lat/lon 76/0, radius 15 degree, point separation 0.5 degree (14.1 KB) Peter, 2017-04-28 10:55

fig.ps (240 KB) Peter, 2017-04-28 11:24

Associated revisions

Revision 18066
Added by Paul 7 months ago

Address issue #1085


#1 Updated by Paul 7 months ago

  • Status changed from New to Feedback

You are missing the -fg option to tell gmtselect that these are not Cartesian y-x coordinates but geographic lat-lon. Adding -fg seems to produce the desired result.

#2 Updated by Peter 7 months ago

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
90.00000 0.00000
$~> echo "90 0" | gmtselect -: -F75_0_15_05.poly -fg
90 0

I guess these are minor issues that may not arise very frequently

#3 Updated by Paul 7 months ago

Thanks, I will look at those. In principle, many programs can read from the same file simultaneously but wonder if there is some limitation on the same program opening multiple streams to the same program.

#4 Updated by Paul 7 months ago

Fixed in r18066. Polygons just touching a pole still have sum(delta_lon) == 0 so a special check was needed. Your cases now pass. I also reworded the introduction to gmtselect to clarify when -fg is need. Thanks for the tests.

#5 Updated by Paul 7 months ago

  • Status changed from Feedback to Resolved
  • Assignee set to Paul
  • Target version set to Candidate for next minor release
  • % Done changed from 0 to 100

Forgot to update settings...

#6 Updated by Joaquim 6 months ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF