[psxy] Plotting a polygon that surrounds the pole

Hello GMT gurus,

I’m running GMT on MacOS 10.14:

gmt psxy [core] 6.1.0_8580087_2020.02.15 [64-bit] - Plot lines, polygons, and symbols in 2-D

I am attempting to plot a polygon that encompasses the South Pole in an orthographic projection and having a little trouble getting GMT to do this for me. Here’s an extract of the relevant code:

psbasemap -Rg -JG150/10/10 -B+gslategray1 -X2 -Y2 -P -K > $output

psxy $icesheets -R -J -G$iceblue -fg -A -O -P >> $output

Running this script fails with the following error:

psxy(72558,0x109e8b5c0) malloc: *** error for object 0x7fadf8000000: pointer being freed was not allocated
psxy(72558,0x109e8b5c0) malloc: *** set a breakpoint in malloc_error_break to debug
lgmpoly_ortho.gmt: line 55: 72558 Abort trap: 6           psxy $icesheets -R -J -G$iceblue -fg -A -O -P -K >> $output

The offending GMT multisegment polygon file, “$icesheets” in my script, can be downloaded here (15 Mb).

Thinking that the problem could be with having a polygon that covers the pole, I also tried splitting the polygon up into 10x10 degree blocks. The plotting still fails, but now with a different error:

Assertion failed: ((fabs(A) > 5 * DBL_EPSILON) || (fabs(B) > 5 * DBL_EPSILON)), function doubleAlmostEqualUlps, file /Users/jkaplan/Downloads/gmt/gmt/src/common_math.c, line 208.
lgmpoly_ortho.gmt: line 56: 72593 Abort trap: 6           psxy $icesheets -R -J -G$iceblue -fg -A -O -P -K >> $output

Any ideas you might have to solve this problem would be much appreciated
Thanks in advance for your help!

Jed

If you could make avbailable some bits that lets me reproduce the problem then happy to give it a go. GMT usually plots polar caps well. Are they true polar cap or are their artificial excursions to the pole one or more points having coordinates equal to that of the pole?

Thanks for the reply Paul,

A link to my data were in the original post, but here it is again, just in case:

Thanks again,

Jed

Sorry, missed that link. I hope to have a look over the weekend, but it is pretty booked already…

No rush. Thanks Paul!

Just a quick gmt info -fg shows there are S poles in there, and a chunk of the data looks like this:

180.0 -84.916702270507812 0
180 -85 0
180.0 -85.003997802734375 0
180.0 -85.087303161621094 0
180.0 -85.170600891113281 0
180.0 -85.253898620605469 0
180.0 -85.337196350097656 0
180.0 -85.420501708984375 0
180.0 -85.503799438476562 0
180.0 -85.58709716796875 0
180.0 -85.670402526855469 0
180.0 -85.753700256347656 0
180.0 -85.836997985839844 0
180.0 -85.920303344726562 0
180.0 -86.00360107421875 0
180.0 -86.086898803710938 0
180.0 -86.170196533203125 0
180.0 -86.253501892089844 0
180.0 -86.336799621582031 0
180.0 -86.420097351074219 0
180.0 -86.503402709960938 0
180.0 -86.586700439453125 0
180.0 -86.669998168945312 0
180.0 -86.753303527832031 0
180.0 -86.836601257324219 0
180.0 -86.919898986816406 0
180.0 -87.003196716308594 0
180.0 -87.086502075195312 0
180.0 -87.1697998046875 0
180.0 -87.253097534179688 0
180.0 -87.336402893066406 0
180.0 -87.419700622558594 0
180.0 -87.502998352050781 0
180.0 -87.5863037109375 0
180.0 -87.669601440429688 0
180.0 -87.752899169921875 0
180.0 -87.836196899414062 0
180.0 -87.919502258300781 0
180.0 -88.002799987792969 0
180.0 -88.086097717285156 0
180.0 -88.169403076171875 0
180.0 -88.252700805664062 0
180.0 -88.33599853515625 0
180.0 -88.419296264648438 0
180.0 -88.502601623535156 0
180.0 -88.585899353027344 0
180.0 -88.669197082519531 0
180.0 -88.75250244140625 0
180.0 -88.835800170898438 0
180.0 -88.919097900390625 0
180.0 -89.002403259277344 0
180.0 -89.085700988769531 0
180.0 -89.168998718261719 0
180.0 -89.252296447753906 0
180.0 -89.335601806640625 0
180.0 -89.418899536132812 0
180.0 -89.502197265625 0
180.0 -89.585502624511719 0
180.0 -89.668800354003906 0
180.0 -89.752098083496094 0
180.0 -89.835403442382812 0
180.0 -89.918701171875 0
180 -90 0
179.998001098632812 -90.0 0
179.914703369140625 -90.0 0
179.831405639648438 -90.0 0
179.748092651367188 -90.0 0
179.664794921875 -90.0 0

which does not look like data but fake lines connecting the perimeter to the S pole. It may be that some tools like arcGIS or Google Earth requires thes fake line to function but for GMT they are not needed and just ends up being taken literally. You could try to search for these sections (probably just one) and remove the junk.

I was having a similar problem and sloved it with a simple awk command before plotting.

#remove the junk S poles
awk ‘$2 > -89.9’ infile.xy > infile_sPole_rm.xy

#plot
psxy nfile_sPole_rm.xy -R -J -Ggrey -fg -A -O -P >> output.ps