gmt
gmt copied to clipboard
WIP Improve detection of closed but split contours
Description of proposed changes
See #6080 for background. After much forensics I learned that the failing contour reported a clockwise (CW) polygon whereas all others are counter-clockwise (CCW). Due to how the contour-trace algorithm and later clipping and projection works, I believe they are all counter-clockwise or should be CCW. Adding a sanity check to converts CW to CCW and prints a warning was added to both grdcontour and pscontour. The correction fixes the problem reported in #6080 and all tests pass in master, except grdcontour/closed.sh gave a slightly different plot (as the change affected what is considered beginning and end. Instead of updating the PS I just upped the RMS limit to 0.02 so the test passes.
Here are all the changes in this PR:
- Change the internal named constants for CCW and CW to be -1 and +1 which better match the uses of handedness in GMT.
- Modify a place in gmt_plot.c where an unset handedness was indicated by GMT_NOTSET. Now we must use 0 instead.
- Add the check, message, and switch if CW is found.
While I am pretty sure this fix is robust (the combined tests found something like ~80 CCW situations for closed contours that were to be ticked and just 1 (in grdcontour/closed.sh) plus the case in #6080 were CW, I am flagging this PR as WIP until we run some more contours with -T+a on more grids.
Not out of the woods yet. Feel free to try this script and run it repeatedly for random locations of the center of the grid:
gmt begin grdcontour png
lon=$(gmt math -Q -180 180 RAND =)
lat=$(gmt math -Q -80 80 RAND =)
gmt grdmath -Rg -I1 $lon $lat SDIST = tmp_dist.nc
gmt coast -Rd -JW$lon/10c -G200 -Sdodgerblue2 -B+t"($lon, $lat)"
gmt grdcontour tmp_dist.nc -A1000+f6p -Wthinnest -T+a -GLZ-/$lon/0,$lon/0/Z+
gmt end show
In general, it looks better but once in a while we have contours (very close to a pole?) and it comes out the wrong way.
I added annotations so we can report bad cases and know which contour it was. Some observations so far:
- Some contours that enclose small areas and wrapped on the E/W sides gets a crazy number of ticks on the East part of the polygon. I think it is always on the east (or xmax). My theory is that somehow the length from the west to east (not drawn) are included in the length estimate which is sued to determined how many ticks.
- We find consistently many contours that are clockwise and I am forcing them to be counter-clockwise (this was the above fix).
- The failing contours are not "closed" on the map (entry and exit at different sides) while the ones that form a closed loop on the map are correctly ticked.
Hi @Esteban82 and @meghanrjones please post if you find other systematic issues.
I have just got this. Is it ok? I don't like the tick in the bottom left and in the upper right just in the borders.
And here the 5000 km contour has ticks in the opposite side.
That one is OK but if you run the script with random lon/lat you will quickly get cases that still fail and give the wrong direction.