gmt icon indicating copy to clipboard operation
gmt copied to clipboard

WIP Improve detection of closed but split contours

Open PaulWessel opened this issue 3 years ago • 6 comments

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:

  1. Change the internal named constants for CCW and CW to be -1 and +1 which better match the uses of handedness in GMT.
  2. Modify a place in gmt_plot.c where an unset handedness was indicated by GMT_NOTSET. Now we must use 0 instead.
  3. 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.

PaulWessel avatar Dec 03 '21 02:12 PaulWessel

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.

PaulWessel avatar Dec 09 '21 22:12 PaulWessel

I added annotations so we can report bad cases and know which contour it was. Some observations so far:

  1. 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.
  2. We find consistently many contours that are clockwise and I am forcing them to be counter-clockwise (this was the above fix).
  3. 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.

PaulWessel avatar Dec 09 '21 22:12 PaulWessel

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.

grdcontour_1

And here the 5000 km contour has ticks in the opposite side.

grdcontour

Esteban82 avatar Dec 12 '21 03:12 Esteban82

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.

PaulWessel avatar Dec 12 '21 04:12 PaulWessel