Bug #851

gmtspatial -C -R clipping

Added by John over 1 year ago. Updated over 1 year ago.

Status:ClosedStart date:2016-01-29
Priority:NormalDue date:
Assignee:Paul% Done:

100%

Category:-
Target version:Candidate for next bugfix release
Affected version:5.x-svn Platform:Mac OS X

Description

Just wanted to officially report a bug in gmtspatial, even though a report was posted to the GMT Forum.

Please test the following...

Given a four-sided polygon (tpoly.xy):

-300 -3500
-200 -800
400 780
500 -3400
-300 -3500

When the command is issued to clip using a -R parameter that is completely inside tpoly.xy:

gmtspatial tpoly.xy -C -R0/100/-3100/-3000

The result comes back with:

21.4815 -3100
27.504 -3000
100 -3000
100 -3100
21.4815 -3100

Expected values are:

0 -3100
0 -3000
100 -3000
100 -3100
0 -3100

Something in the computation for the x-coordinates of the left side of the resulting polygon is in error.

Thanks for resolving.

Associated revisions

Revision 15624
Added by Paul over 1 year ago

Address issue #851

History

#1 Updated by John over 1 year ago

The tpoly.xy should be (the 780 Y-coord should have been negative):

-300 -3500
-200 -800
400 -780
500 -3400
-300 -3500

For the example, above, it probably doesn't matter. Using -780 in tpoly.xy, gmtspatial tpoly.xy -C -R0/100/-3100/-3000 still gives an erroneous result:

21.4815 -3100
27.504 -3000
100 -3000
100 -3100
21.4815 -3100

Using the same tpoly.xy, but changing the -R parameters to be negative in X (e.g., -100/0/-3100/-3000) yields a null result (i.e., just sends a line feed).

Changing -R to 150/250/-3100/-3000) yields a correct result. Makes me wonder whether things are being considered as geographic lon/lats, when they're actually Cartesian?

#2 Updated by John over 1 year ago

New example. Using GMT5-svn version r15528:

First, let tpoly.xy be a square of size 100, with lower left corner at the origin:

0 0
0 100
100 100
100 0
0 0

gmtspatial tpoly.xy -C -R0/10/0/10 (lower left corner) yields an proper result (two shared sides):

0 0
10 0
10 10
0 10
0 0

Now, increase tpoly.xy to be a square of size 200:

0 0
200 0
200 200
0 200
0 0

and run the same command: gmtspatial tpoly.xy -C -R0/10/0/10, an erroneous result is given:

0 0
0 0
0 10
0 10
0 0

The breaking point occurs when the size of the polygon square is increased past 190. Square tpoly.xy of 190 "works;" square tpoly.xy of 191, "fails." Very odd, indeed...

The fact that 190 is 180 past the eastern limit of 10, makes me think that something in the business of "date line" is messing things up.

#3 Updated by John over 1 year ago

Line 1574 in gmtspatial has:

for (p = 0; p < np; p++) GMT_xy_to_geo (GMT, &cp[GMT_X][p], &cp[GMT_Y][p], cp[GMT_X][p], cp[GMT_Y][p]);

Is this being applied to non-geographic values, too? Culprit?

#4 Updated by Paul over 1 year ago

  • Status changed from New to In Progress

I have added your two failing examples as new test scripts for gmtspatial under 5.2.2 (and trunk). They both fail but we will now work on fixing the issues.

#5 Updated by John over 1 year ago

Thanks, Paul! I will be happy to test it with some "real world" applications, once you have an update. My own work-around has been a brute force approximation using high-resolution grids. When there is a fix, my software will be a lot cleaner and the answers should be more rigorous. I'll keep an eye on an e-mail letting me know the progress. I know you've been busy! Thanks, again!

#6 Updated by Paul over 1 year ago

  • Status changed from In Progress to Resolved
  • Assignee set to Paul
  • Target version set to Candidate for next bugfix release
  • % Done changed from 0 to 100

Problem was an insistence on calling GMT_wesn_clip instead of GMT_rect_clip when data are Cartesian... Fixed in r15626.

#7 Updated by John over 1 year ago

After svn update to r15639, I can confirm that gmtspatial works as expected in real world examples.

Thanks again, Paul!

#8 Updated by Joaquim over 1 year ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF