Bug #497

Projection error when minimum longitude is -180

Added by Eduardo almost 4 years ago. Updated almost 4 years ago.

Status:ClosedStart date:2014-01-31
Priority:NormalDue date:
Assignee:-% Done:

0%

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

Description

Data is -180/180 in longitude. I've tested cylindrical equidistant (with lat=43.5) and Winkel Tripel.

Cylindrical Equidistant


psxy land_polygons.gmt -R-180/-28/-56/84 -JQ-103/43.5/25h -Ggreen@50 -W0.1p,black,solid -P > america2.ps
ERROR: Caught signal number 11 (Segmentation fault) at
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_cyleqdist+0x52)[0x7f94c1961892]
[0x7f94aa8b2000]
Stack backtrace:
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(sig_handler+0x13a)[0x7f94c191580a]
/lib/x86_64-linux-gnu/libpthread.so.0(+0xf210)[0x7f94bfa79210]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_cyleqdist+0x52)[0x7f94c1961892]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_geo_to_xy+0x74)[0x7f94c19e3b34]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_move_to_wesn+0x83)[0x7f94c19e4de3]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_wesn_clip_old+0x25a)[0x7f94c19e8f3a]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_wesn_clip+0x1059)[0x7f94c19ee989]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_clip_to_map+0xea)[0x7f94c19e86ca]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_geo_polygon+0x95)[0x7f94c1a0fca5]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_geo_polygon_segment+0x34)[0x7f94c1a100f4]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_geo_polygons+0x105)[0x7f94c1a10345]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_psxy+0x22af)[0x7f94c1b1cd9f]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_Call_Module+0xab)[0x7f94c19253ab]
psxy(main+0x27d)[0x400f7d]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5)[0x7f94bf6df995]
psxy[0x400c39]

Winkel Tripel


psxy land_polygons.gmt -R-180/-28/-56/84 -JR-103/25h -Ggreen@50 -W0.1p,black,solid -P > america.ps
ERROR: Caught signal number 11 (Segmentation fault) at
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_winkel+0x121)[0x7fb6cf79e341]
[0x7fb6b90e8000]
Stack backtrace:
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(sig_handler+0x13a)[0x7fb6cf74880a]
/lib/x86_64-linux-gnu/libpthread.so.0(+0xf210)[0x7fb6cdaba210]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_winkel+0x121)[0x7fb6cf79e341]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_geo_to_xy+0x74)[0x7fb6cf816b34]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_move_to_wesn+0x83)[0x7fb6cf817de3]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_wesn_clip_old+0x25a)[0x7fb6cf81bf3a]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_wesn_clip+0x1059)[0x7fb6cf821989]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_clip_to_map+0xea)[0x7fb6cf81b6ca]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_geo_polygon+0x95)[0x7fb6cf842ca5]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(gmt_geo_polygon_segment+0x34)[0x7fb6cf8430f4]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_geo_polygons+0x105)[0x7fb6cf843345]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_psxy+0x22af)[0x7fb6cf94fd9f]
/home/opt/apps/GMT/GMT-5-testing/bin/../lib/x86_64-linux-gnu/libgmt.so.5(GMT_Call_Module+0xab)[0x7fb6cf7583ab]
psxy(main+0x27d)[0x400f7d]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5)[0x7fb6cd720995]
psxy[0x400c39]

05-02-2014 16-13-44.jpg (160 KB) Joaquim, 2014-02-05 06:19

to_paul.bz2 (546 KB) Eduardo, 2014-02-05 16:35

eurasia.png (21.1 KB) Joaquim, 2014-02-06 04:56

History

#1 Updated by Joaquim almost 4 years ago

Eduaro,
On Windows when I do the equivalent of

echo 0 | gawk '{printf "-180 0\n0 0\n180 1"}' | psxy -R-180/-28/-56/84 -JQ-103/43.5/25h > lixo.ps

it works well. Does it work for you too?

#2 Updated by Eduardo almost 4 years ago

Yes, that works for me too.

My datafile is the entire Openstreetmap land polygons. N = 32031744 data points. The file size is 680M. My machine has 16GB RAM.

With this datafile it works using -R-179/... instead of R-180/...

If I do a gmtselect before running psxy it also works.

Here are a couple of backtraces from gdb. Both have problem at the same data point, when j=2097150

Cylindrical Equidistant

#0  0x00007ffff7a0b249 in GMT_cyleqdist (GMT=0x7fffe7b10010, lon=75, 
    lat=44.964836599999998, x=0x7fffe092f000, y=0x7fffdf92e000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:283
#1  0x00007ffff7ae3669 in GMT_geo_to_xy (GMT=0x7fffe7b10010, lon=-28, 
    lat=44.964836599999998, x=0x7fffe092f000, y=0x7fffdf92e000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:2604
#2  0x00007ffff7adc804 in gmt_move_to_wesn (GMT=0x7fffe7b10010, 
    x_edge=0x7fffdf92f010, y_edge=0x7fffde92e010, lon=-28, 
    lat=44.964836599999998, lon_old=37.277333599999999, 
    lat_old=44.962180699999998, j=2097150, nx=0)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1472
#3  0x00007ffff7adcb99 in gmt_wesn_clip_old (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, n=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8, total_nx=0x7fffffff1320)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1503
#4  0x00007ffff7adce3c in GMT_wesn_clip (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, n_orig=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8, total_nx=0x7fffffff1320)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1556
#5  0x00007ffff7ada014 in GMT_clip_to_map (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, np=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1100
#6  0x00007ffff7b25d10 in gmt_geo_polygon (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, n=4358693)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4090
#7  0x00007ffff7b26528 in gmt_geo_polygon_segment (GMT=0x7fffe7b10010, 
    PSL=0x40bda0, S=0xb6e4450, add_pole=false)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4186
#8  0x00007ffff7b267dc in GMT_geo_polygons (GMT=0x7fffe7b10010, S=0xb6e4450)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4236
#9  0x00007ffff7cee20e in GMT_psxy (V_API=0x40ae90, mode=6, 
    args=0x7fffffffe760)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/psxy.c:1144
#10 0x00007ffff79b08c6 in GMT_Call_Module (V_API=0x40ae90, 
    module=0x412710 "psxy", mode=6, args=0x7fffffffe760)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_api.c:5893
#11 0x0000000000401587 in main (argc=7, argv=0x7fffffffe758)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt.c:178

Winkel Tripel

#0  0x00007ffff7a1baa8 in GMT_winkel (GMT=0x7fffe7b10010, 
    lon=0.6544984694978736, lat=0.78478444629125244, x=0x7fffe092f000, 
    y=0x7fffdf92e000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:2096
#1  0x00007ffff7ae3669 in GMT_geo_to_xy (GMT=0x7fffe7b10010, lon=-28, 
    lat=44.964836599999998, x=0x7fffe092f000, y=0x7fffdf92e000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:2604
#2  0x00007ffff7adc804 in gmt_move_to_wesn (GMT=0x7fffe7b10010, 
    x_edge=0x7fffdf92f010, y_edge=0x7fffde92e010, lon=-28, 
    lat=44.964836599999998, lon_old=37.277333599999999, 
    lat_old=44.962180699999998, j=2097150, nx=0)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1472
#3  0x00007ffff7adcb99 in gmt_wesn_clip_old (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, n=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8, total_nx=0x7fffffff1320)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1503
#4  0x00007ffff7adce3c in GMT_wesn_clip (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, n_orig=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8, total_nx=0x7fffffff1320)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1556
#5  0x00007ffff7ada014 in GMT_clip_to_map (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, np=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1100
#6  0x00007ffff7b25d10 in gmt_geo_polygon (GMT=0x7fffe7b10010, 
    lon=0x7fffccffc010, lat=0x7fffd0ffd010, n=4358693)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4090
#7  0x00007ffff7b26528 in gmt_geo_polygon_segment (GMT=0x7fffe7b10010, 
    PSL=0x40bda0, S=0xb6e4450, add_pole=false)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4186
#8  0x00007ffff7b267dc in GMT_geo_polygons (GMT=0x7fffe7b10010, S=0xb6e4450)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4236
#9  0x00007ffff7cee20e in GMT_psxy (V_API=0x40ae90, mode=6, 
    args=0x7fffffffe760)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/psxy.c:1144
#10 0x00007ffff79b08c6 in GMT_Call_Module (V_API=0x40ae90, 
    module=0x412710 "psxy", mode=6, args=0x7fffffffe760)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_api.c:5893
#11 0x0000000000401587 in main (argc=7, argv=0x7fffffffe758)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt.c:178

Here is an extract of the data file, nothing strange.

Record #    longitude     latitude
2097140    42.8413677    12.9668687
2097141    42.8416325    12.9666133
2097142    42.8421588    12.9664894
2097143    42.8422985    12.9654603
2097144    42.8424348    12.9649464
2097145    42.8421709    12.9650724
2097146    42.841382     12.9651951
2097147    42.8413843    12.9649375
2097148    42.8411229    12.9648061
2097149    42.8409921    12.9646767
2097150    42.8407341    12.9641596
2097151    42.8406043    12.9640291
2097152    42.8404737    12.9638999
2097153    42.840344     12.9637695
2097154    42.8398211    12.9635075
2097155    42.8395584    12.9635054
2097156    42.9410587    12.9795245
2097157    42.9405314    12.979778
2097158    42.9403986    12.9799061
2097159    42.9402666    12.9800335
2097160    42.9389523    12.9801522

As I said before, it works with a data subset, even with a large one (N = 2099855).

#3 Updated by Paul almost 4 years ago

  • Status changed from New to Feedback

Hi Eduardo-

I made the ogr output and running psxy and it has been running for a while now; will let you know if it ever finishes or crashes. This is OSX though so I bet it will eventually work...

#4 Updated by Paul almost 4 years ago

I added some -Vl feedback to psxy (in r12819). What I find is that for me, psxy reads the data (takes a minute or so), then races through the segments until it hits
psxy: Plotting table 0 segment 32391
and there it sits forever. Please let me know if this is where you crash (update GMT first).
-p

#5 Updated by Paul almost 4 years ago

It hangs on the Eurasia polygon which goes from lon = -17.5301717 to 180 because if is truncated at lon = 180. Presumably the far Siberia bit is a separate polygon with lonts starting at -180. When I killed my psxy it says it was in fixup_path, probably to do with that part. I am rerunning with -A; perhaps you can try that too.
-p

#6 Updated by Paul almost 4 years ago

WIth -A the fixup_path is bypassed but now spending lots of time in psl_shorten_path. I am just plotting the EurAfricasia polygon and still running in debug...

#7 Updated by Eduardo almost 4 years ago

psxy land_polygons.gmt -R-180/-28/-56/84 -JR-103/25h -Ggreen@50 -W0.1p,black,solid -P -Vl > america.ps

...
psxy: Plotting table 0 segment 209500
psxy: Plotting table 0 segment 209501
psxy: Plotting table 0 segment 209502
ERROR: Caught signal number 11 (Segmentation fault)...
psxy land_polygons.gmt -R-20/195/10/84 -JR90/20 -Ggreen@50 -W0.1p,black,solid -Vl > eurasia.p

...
psxy: Plotting table 0 segment 451138
psxy: Plotting table 0 segment 451139
psxy: Plotting table 0 segment 451140
psxy: Plotting table 0 segment 451141
psxy: Plotting table 0 segment 451142
psxy: Plotting table 0 segment 451143
^C

no segfault, break at psl_shorten_path

There is a small hiccup at segment 209502 and then continues.

I did a strace of the running process, and before the stall it says the same that at 209502:

write(2, "psxy: Plotting table 0 segment 4"..., 38) = 38
mmap(NULL, 33558528, PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANONYMOUS, -1, 0) = 0x7fe6060e8000
mmap(NULL, 33558528, PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANONYMOUS, -1, 0) = 0x7fe6040e7000
munmap(0x7fe5fd13d000, 28876800)        = 0
munmap(0x7fe60ad86000, 28876800)        = 0
mremap(0x7fe6060e8000, 33558528, 28876800, MREMAP_MAYMOVE) = 0x7fe6060e8000
mremap(0x7fe6040e7000, 33558528, 28876800, MREMAP_MAYMOVE) = 0x7fe6040e7000
mmap(NULL, 33558528, PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANONYMOUS, -1, 0) = 0x7fe60ad88000
mmap(NULL, 33558528, PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANONYMOUS, -1, 0) = 0x7fe5fdd0c000
mremap(0x7fe60ad88000, 33558528, 28876800, MREMAP_MAYMOVE) = 0x7fe60ad88000
mremap(0x7fe5fdd0c000, 33558528, 28876800, MREMAP_MAYMOVE) = 0x7fe5fdd0c000

with -A I get a nice segfault at the same line

psxy land_polygons.gmt -R-20/195/10/84 -JR90/20 -A -Ggreen@50 -W0.1p,black,solid -Vl > eurasia.ps

...
psxy: Plotting table 0 segment 451138
psxy: Plotting table 0 segment 451139
psxy: Plotting table 0 segment 451140
psxy: Plotting table 0 segment 451141
psxy: Plotting table 0 segment 451142
psxy: Plotting table 0 segment 451143
ERROR: Caught signal number 11 (Segmentation fault) at...

#8 Updated by Paul almost 4 years ago

and if you extract segment 451143 with gmtconvert -Q451143 land_polygons.gmt > junk.txt and run that by itself, do you get the same error?

#9 Updated by Eduardo almost 4 years ago

yes!

psxy junk.txt -R-20/195/10/84 -JR90/20 -A -Ggreen@50 -W0.1p,black,solid -Vl > eurasia.pspsxy: Processing input table data
psxy: Warning: spherical approximation used!
psxy: Reading Data Table from File junk.txt
psxy: Plotting table 0 segment 0
ERROR: Caught signal number 11 (Segmentation fault)

This is the continental Americas polygon:

junk.txt: N = 3609105    <-168.1239999/-34.7930182>    <-53.89706/72.0015>

The weird thing is that I can plot it with the right -R.

#10 Updated by Paul almost 4 years ago

The size of these polygons makes debug painful - I am still running the Eurafricasia polygon.
Perhaps try the junk.txt after decimating, e.g., first only use every other point, if fail, try every 4th point etc until you have the smallest file that fails. Might be easier then to load in debugger and see what exactly is happening.

#11 Updated by Eduardo almost 4 years ago

well, using a junk.txt subset the error is at point

-64.9033944 46.8991115

The 2097096th point. (The file has 2097097 points, the first is >.)

If I delete one data point randomly (not the >) in the same file, psxy runs OK.

#12 Updated by Paul almost 4 years ago

Now that is silly. And cannot tell what line it crashes at?

#13 Updated by Paul almost 4 years ago

I mean what line in GMT it crashes at?

#14 Updated by Joaquim almost 4 years ago

Since MS debugger is a superior tool, can you make the file available to me so that I can try my (our) luck?

#15 Updated by Eduardo almost 4 years ago

jajaja

it's a 212 MB file.

http://fcaglp.fcaglp.unlp.edu.ar/~esuarez/issue_497/land_polygons.gmt.bz2

to me is related to some memory allocation problem.

I've narrowed it to this scenario:

  • You have a HUGE data file.
  • You want to paint the polygons, drawing only the lines doesn't trigger
    the problem.
  • the plot region includes a tiny part of another HUGE polygon. In my case, the
    problem arises in the Bering Sea because I'm ploting Eurasia or the Americas
    so the plot includes a tiny part of the other continent there.
    An experiment I did was try to plot Alaska using the entire file or gmtselecting
    only the region of interest. It failed with the entire file but not with the
    sub-region.

#16 Updated by Paul almost 4 years ago

Does adding -A avoid the crash?

#17 Updated by Eduardo almost 4 years ago

No, using -A does anything.

As an example, to me this fails:

psxy land_polygons.gmt -R-190/-60/50/90 -JQ-125.5/43.5/15 -Glightgreen -W0.25p,black,solid -P -Vl > test_1.ps

And this works:

gmtselect land_polygons.gmt -R-190/-60/50/90 > north.d
psxy north.d -R-190/-60/50/90 -JQ-125/43.5/15 -Glightgreen -W0.25p,black,solid -P -Vl > test_2.ps

#18 Updated by Joaquim almost 4 years ago

I might be getting close. At line gmt_map line 1103, x & y are undetermined (never got initialized) and a crash ends up happening in GMT_winkel() when accessing to this pointers.

else    /* Mixed case so we must clip the polygon */
    n = (*GMT->current.map.clip) (GMT, lon, lat, np, x, y, &total_nx);

#19 Updated by Paul almost 4 years ago

Hm, okey. so what is np etc and from where is this called (i.e., there should not be possible to pass down x,y pointer-pointers that did not point to actual memory)

#20 Updated by Joaquim almost 4 years ago

out = 3596752
np = 3609105

#21 Updated by Eduardo almost 4 years ago

Joaquim, the problem isn't exclusive of Winkel. Whereas Cylindrical and Miscellaneous fail, Conic and azimuthal work.

Here are where the routines gives the segmentation fault. They fail at the same instruction, when they try to asign something to *x

Program received signal SIGSEGV, Segmentation fault.
0x00007ffff7a0b904 in GMT_cylstereo (GMT=0x7fffe810a010, lon=65.5, lat=50, 
    x=0x7fffe0f29000, y=0x7fffdff28000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:339
339             *x = lon * GMT->current.proj.j_x;

Program received signal SIGSEGV, Segmentation fault.
0x00007ffff7a0adb0 in GMT_cyleq (GMT=0x7fffe810a010, lon=65.5, 
    lat=49.873610229086331, x=0x7fffe0f29000, y=0x7fffdff28000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:246
246             *x = lon * GMT->current.proj.j_x;

Program received signal SIGSEGV, Segmentation fault.
0x00007ffff7a1ee77 in GMT_sinusoidal (GMT=0x7fffe810a010, lon=65.5, 
    lat=0.87045870835387995, x=0x7fffe0f29000, y=0x7fffdff28000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:2496
2496            *x = GMT->current.proj.EQ_RAD * lon * D2R * cos (lat);

Program received signal SIGSEGV, Segmentation fault.
0x00007ffff7a1af0b in GMT_grinten (GMT=0x7fffe810a010, lon=65.5, lat=50, 
    x=0x7fffe0f29000, y=0x7fffdff28000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:2028
2028            *x = copysign (GMT->current.proj.v_r, lon) * (A * GP2 + sqrt (A2 * GP2 * GP2 - P2A2 * (G*G - P2))) * i_P2A2;

Program received signal SIGSEGV, Segmentation fault.
0x00007ffff7a1e7ab in GMT_robinson (GMT=0x7fffe810a010, lon=65.5, lat=50, 
    x=0x7fffe0f29000, y=0x7fffdff28000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:2427
2427            *x = GMT->current.proj.n_cx * X * lon;  /* D2R is in n_cx already */

#22 Updated by Joaquim almost 4 years ago

Eduardo, yes I know is not a Winkel specific. The problem is what I identified in comment #18. The x and y variables are in that case pointing to NULL and naturally a crash follows.

#23 Updated by Paul almost 4 years ago

Can you tell me from where this is called when x,y are NULL? I.e., we start at GMT_geo_polygons I assume and then there is perhaps clipping and wrapping and we end up with a piece that somehow has NULL?

#24 Updated by Joaquim almost 4 years ago

Yes, we start a gmt_geo_polygon line 4091 and than at GMT_clip_to_map if out is neither equal zero or equal np the third path always use uninitialized x & y

#25 Updated by Paul almost 4 years ago

Not so, I think. In gmt_geo_polygon, you have double *xp = NULL, *yp = NULL; so both are initialized as NULL pointers. Then their ADDRESSES (which cannot be NULL) are passed to GMT_clip_to_map so that they can return pointers to memory allocated deeper down. Thus inside GMT_clip_to_map the x,y in the last case will always have some address. The question is where the memory allocation fails. I assume the pointer to the clip routine takes you to GMT_wesn_clip, and inside there the *x and *y are allocated. I suspect there is some case where not enough memory is allocated and you try to write to that segment. For me under OSX it does not crash so hard to check.

#26 Updated by Paul almost 4 years ago

Just to clarify: Since Eduardo can plot it if his -R does not induce clipping, it is clearly the map clipping that is failing for him. So when things go boom, can you tell what the calling stack is? I understand it crashes at some line like

*x = GMT→current.proj.n_cx * X * lon;

but what I need to know is which function call that, and which one called that, etc. so I can get a picture of where in the algorithm you are. Can you step up and see if that address being passed to the lowest gmt_proj function is out of range?

#27 Updated by Paul almost 4 years ago

Rereading Eduardo's posting #2 I examined gmt_wesn_clip_old and found an issue that might be the problem. In r12854; let me know if it made any difference at all.

#28 Updated by Eduardo almost 4 years ago

Fail in GMT 5.1.1 r12855. Here is the backtrace from gdb:

#0  0x00007ffff7a1bb50 in GMT_winkel (GMT=0x7fffe7ec6010, 
    lon=0.6544984694978736, lat=0.78478444629125244, x=0x7fffe0ce5000, 
    y=0x7fffdfce4000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_proj.c:2096
#1  0x00007ffff7ae38a7 in GMT_geo_to_xy (GMT=0x7fffe7ec6010, lon=-28, 
    lat=44.964836599999998, x=0x7fffe0ce5000, y=0x7fffdfce4000)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:2617
#2  0x00007ffff7adca42 in gmt_move_to_wesn (GMT=0x7fffe7ec6010, 
    x_edge=0x7fffdfce5010, y_edge=0x7fffdece4010, lon=-28, 
    lat=44.964836599999998, lon_old=37.277333599999999, 
    lat_old=44.962180699999998, j=2097150, nx=0)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1485
#3  0x00007ffff7adcdd7 in gmt_wesn_clip_old (GMT=0x7fffe7ec6010, 
    lon=0x7fffcd3b2010, lat=0x7fffd13b3010, n=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8, total_nx=0x7fffffff1320)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1516
#4  0x00007ffff7add07a in GMT_wesn_clip (GMT=0x7fffe7ec6010, 
    lon=0x7fffcd3b2010, lat=0x7fffd13b3010, n_orig=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8, total_nx=0x7fffffff1320)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1569
#5  0x00007ffff7ada252 in GMT_clip_to_map (GMT=0x7fffe7ec6010, 
    lon=0x7fffcd3b2010, lat=0x7fffd13b3010, np=4358693, x=0x7fffffff13d0, 
    y=0x7fffffff13c8)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_map.c:1113
#6  0x00007ffff7b25d7c in gmt_geo_polygon (GMT=0x7fffe7ec6010, 
    lon=0x7fffcd3b2010, lat=0x7fffd13b3010, n=4358693)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4089
#7  0x00007ffff7b26594 in gmt_geo_polygon_segment (GMT=0x7fffe7ec6010, 
    PSL=0x40bd80, S=0xb6e4600, add_pole=false)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4185
#8  0x00007ffff7b26848 in GMT_geo_polygons (GMT=0x7fffe7ec6010, S=0xb6e4600)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_plot.c:4235
#9  0x00007ffff7cee29c in GMT_psxy (V_API=0x40ae70, mode=6, 
    args=0x7fffffffe760)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/psxy.c:1150
#10 0x00007ffff79b096f in GMT_Call_Module (V_API=0x40ae70, 
    module=0x412850 "psxy", mode=6, args=0x7fffffffe760)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt_api.c:5978
#11 0x0000000000401597 in main (argc=7, argv=0x7fffffffe758)
    at /root/software/gmt/GMTdev/gmt5/trunk/src/gmt.c:178

#29 Updated by Joaquim almost 4 years ago

I did navigate back and forth along the stack trace but cannot find an obvious place where that allocation would take place. If I had a single segment causing the crash I could do a step by step execution around the trouble zone but as is I have to wait till it crash and than investigate.
The best place for the allocation seam to be in gmt_wesn_clip_old() (see attached image) but still cannot understand why memory is invalid after the call at lines 1489 and 1490 but the content of GMT→hidden.mem_coord[GMT_X] is not.

#30 Updated by Paul almost 4 years ago

OK, pretty sure there are issues in gmt_move_to_wesn where we add to an array that apparently is not allocated or set. I will add more debug verbosity in there but busy the next several hours... If you want to try it has to do with how memory is allocated in gmt_wesn_clip_old via GMT_prep_tmp_arrays. Perhaps add some checking there or stop there in debug.

#31 Updated by Joaquim almost 4 years ago

I'm close to find it (or not, because I thought that before) but more eyes looking the better and I'll have to interrupt my debug session soon.
The problem is that GMT→hidden.mem_rows is reset back to GMT_INITIAL_MEM_ROW_ALLOC = 2097152 and since that is not enough to hold the polygon it naturally kabooms. Now the mystery to solve is how is it set back to that value since I followed it in gmt_memory it was correctly allocated to the double of that value.

But since we are here, I see that we are allocating by trying on double. We had discussed this time ago and agreed that this is overkill. We don't want to try to go from 500 Mb to 1 Gb. The double strategy can only go up to a (currently) reasonable value, than it has to proceed at a constant step.

while (row >= GMT->hidden.mem_rows) GMT->hidden.mem_rows <<= 1;    /* Double up until enough */

#32 Updated by Joaquim almost 4 years ago

OK VICTORY

In gmt_wesn_clip_old() the call to

GMT_prep_tmp_arrays (GMT, 1, 2);    /* Init or reallocate tmp vectors */

because of the argument "1" forces a reallocation to the initial (and too short) GMT_INITIAL_MEM_ROW_ALLOC.

Paul, I'll let you do the fix as I would be touching my way out of this.

#33 Updated by Paul almost 4 years ago

Note my fix yesterday was to change that from GMT_prep_tmp_arrays (GMT, 0, 2); to GMT_prep_tmp_arrays (GMT, 1, 2); since we needed at least one - then it should expand as needed when more points are added.

#34 Updated by Joaquim almost 4 years ago

Sorry, bad info. The culprit is not the "1" arg but rather GMT_assign_segment() that at line 5946 sets

GMT->hidden.mem_cols = 0;    /* Flag that we need to reallocate new temp arrays for next segment, if any */

which triggers later down in GMT_prep_tmp_arrays the killing realloc

if (GMT->hidden.mem_cols == 0)
    gmt_init_tmp_arrays (GMT, n_cols);    /* First time we get here */

#35 Updated by Paul almost 4 years ago

Not sure I understand. Which line in GMT_init_tmp_arrays causes a crash? BTW, in r12857 I have added more debug print statements so please rerun case with -Vd and see if output is helpful.

#36 Updated by Joaquim almost 4 years ago

I don't need more debug prints. Problem is identified.

in gmt_wesn_clip_old the call to GMT_prep_tmp_arrays (GMT, 1, 2); causes the condition

if (GMT→hidden.mem_cols == 0)
gmt_init_tmp_arrays (GMT, n_cols); /* First time we get here */

to execute. But this call resets the mem allocation (that had previously be done correctly upon file reading) to the default value of GMT_INITIAL_MEM_ROW_ALLOC, which is too small and this is reason for the crash.

#37 Updated by Paul almost 4 years ago

I understand that at the end of reading the very large segment, GMT_assign_segment is called and it simply passes the temp arrays off to the segment (instead of duplicating it), sets the men_coard[] array to NULL and sets men_cols = 0. Fine. Then we need to clip the array and GMT_prep_tmp_arrays (GMT, 1,2) is called. As you correctly say, this calls gmt_init_tmp_arrays (GMT, 2) to be called since mem_cols == 0. This function reallocates the men_coard array (still to length 2), sets men_cols to 2, and allocates column arrays of length GMT→hidden.mem_rows = GMT_INITIAL_MEM_ROW_ALLOC. Presumably all is well so far. Then gmt_wesn_clip_old continues with the loop over n, with n >> GMT_INITIAL_MEM_ROW_ALLOC. Then as crossings are found we add them to the xx,yy arrays but we call GMT_prep_tmp_arrays (GMT, j, 2); first. This should ensure that j >= GMT→hidden.mem_rows then more memory is allocated first. So when xx[j] = xc[k] happens the xx array is long enough. Then after that we may add 0,1, or 2 points via gmt_move_to_wesn, hence we call GMT_prep_tmp_arrays (GMT, j+2, 2); before that step, with the +2 ensuring we would have enough space. As we loop over i→n we end up allocating more space by doubling GMT→hidden.mem_rows.

So I still do not understand how allocating fresh new arrays of length GMT_INITIAL_MEM_ROW_ALLOC end up being too short at some point. Can you help me understand? Perhaps those print statements I added can help me do that?

#38 Updated by Joaquim almost 4 years ago

OK, here are the printings (my file has 2100000 records)

psxy: Processing input table data
psxy: Warning: spherical approximation used!
psxy: Reading Data Table from File lixo1.dat
psxy: Plotting table 0 segment 0

C:\progs_cygw\GMTdev\gmt5\trunk\build>
C:\progs_cygw\GMTdev\gmt5\trunk\build>
C:\progs_cygw\GMTdev\gmt5\trunk\build>psxy -R-20/195/10/84 -JR90/20 -A -Ggreen@50 -W0.1p,black,solid lixo1.dat -Vl > eurasia.ps
psxy: Processing input table data
psxy: Warning: spherical approximation used!
psxy: Reading Data Table from File lixo1.dat
^C
C:\progs_cygw\GMTdev\gmt5\trunk\build>psxy -R-20/195/10/84 -JR90/20 -A -Ggreen@50 -W0.1p,black,solid lixo1.dat -Vd > eurasia.ps
psxy: Processing input table data
psxy: Warning: spherical approximation used!
psxy: Projected values in meters: -9.93762e+006 9.48675e+006 1.20656e+006 9.46513e+006
psxy: GMTAPI_Init_Import: Passed family = Data Table and geometry = Polygon
psxy: Object ID 0 : Registered Data Table File lixo1.dat as an Input resource with geometry Polygon [n_objects = 1]
==> 1 API Objects at end of GMT_Register_IO
-----------------------------------------------------------
K.. ID RESOURCE.... DATA........ FAMILY.... DIR... S O M L
-----------------------------------------------------------
* 0  0            0            0 Data Table  Input 0 Y N 1
-----------------------------------------------------------
psxy: GMTAPI_Init_Import: Added 1 new sources
psxy: GMT_Init_IO: Returned first Input object ID = 0
psxy: GMTAPI_Begin_IO: Input resource access is now enabled [container]
psxy: GMTAPI_Import_Dataset: Passed ID = -1 and mode = 0
psxy: Reading Data Table from File lixo1.dat
psxy: GMT memory: Initialize 2 temporary column arrays, each of length : 2097152
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304
psxy: GMT_assign_segment: Pass 2 large arrays with length = 2100000 off and get new tmp arrays
psxy: GMT_set_seg_polar: N = 2100000 Multiples of 360: 0  Residual: 0 Polygon contains no pole.
psxy: GMT_End_IO: Input resource access is now disabled
==> 1 API Objects at end of GMT_Get_Data
-----------------------------------------------------------
K.. ID RESOURCE.... DATA........ FAMILY.... DIR... S O M L
-----------------------------------------------------------
  0  0            0       2e24b0 Data Table  Input 2 Y N 1
-----------------------------------------------------------
==> 1 API Objects at end of GMT_Read_Data
-----------------------------------------------------------
K.. ID RESOURCE.... DATA........ FAMILY.... DIR... S O M L
-----------------------------------------------------------
  0  0            0       2e24b0 Data Table  Input 2 Y N 1
-----------------------------------------------------------
psxy: Plotting table 0 segment 0
psxy: GMT memory: Initialize 2 temporary column arrays, each of length : 2097152
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304

#39 Updated by Eduardo almost 4 years ago

Here are the last part of my printings

...
psxy: Plotting table 0 segment 209498
psxy: Plotting table 0 segment 209499
psxy: Plotting table 0 segment 209500
psxy: Plotting table 0 segment 209501
psxy: Plotting table 0 segment 209502
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304
psxy: GMT memory: Increase 2 temporary column arrays to new length : 8388608
psxy: GMT memory: Initialize 2 temporary column arrays, each of length : 2097152
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304
ERROR: Caught signal number 11 (Segmentation fault)...

#40 Updated by Paul almost 4 years ago

Joaquim, regarding post 38, it then crashes after successfully allocating to 4194304? Eduardo, in posting 39, so you never got any "GMT_assign_segment: Pass 2 large arrays with length = xxxxxxx off and get new tmp arrays" message?

#41 Updated by Joaquim almost 4 years ago

OK, another illusion of having found where (not why). The debugger told me that the crash follows when in gmt_wesn_clip_old()

GMT_prep_tmp_arrays (GMT, j+2, 2);    /* Init or reallocate tmp vectors */
j += gmt_move_to_wesn (GMT, xx, yy, lon[i], lat[i], lon[i-1], lat[i-1], j, nx);    /* May add 2 points, which explains the j+2 stuff */

j = GMT_INITIAL_MEM_ROW_ALLOC - 2

that is, when new memory is (re)allocated. And in fact I could verify that at the right moment that happens inside GMT_prep_tmp_arrays()

GMT→hidden.mem_coord[col] = GMT_memory (GMT, GMT→hidden.mem_coord[col], GMT→hidden.mem_rows, double);

the xx and yy pointers IMMEDIATELY get corrupted. So it's the reallocation action that kills.
I know this makes very little sense but one thing strikes me. As we can see from the printings, the memory (the large chunks) seam to be allocated twice. Why? What happened in between?

#42 Updated by Eduardo almost 4 years ago

Only two. Full log attached.

...
psxy: GMT_set_seg_polar: N = 281 Multiples of 360: 0  Residual: 0 Polygon contains no pole.
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304
psxy: GMT memory: Increase 2 temporary column arrays to new length : 8388608
psxy: GMT_assign_segment: Pass 2 large arrays with length = 4358693 off and get new tmp arrays
psxy: GMT_set_seg_polar: N = 4358693 Multiples of 360: 0  Residual: -6.29115e-11 Polygon contains no pole.
psxy: GMT memory: Initialize 2 temporary column arrays, each of length : 2097152
psxy: GMT_set_seg_polar: N = 15372 Multiples of 360: 0  Residual: 0 Polygon contains no pole.
...
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304
psxy: GMT_assign_segment: Pass 2 large arrays with length = 3609105 off and get new tmp arrays
psxy: GMT_set_seg_polar: N = 3609105 Multiples of 360: 0  Residual: 0 Polygon contains no pole.
psxy: GMT memory: Initialize 2 temporary column arrays, each of length : 2097152
psxy: GMT_set_seg_polar: N = 19133 Multiples of 360: 0  Residual: 0 Polygon contains no pole
...

#43 Updated by Paul almost 4 years ago

For posting #41:

When GMT starts, the temp arrays are initially allocated to 2097152, always. then

psxy: Plotting table 0 segment 209502
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304 [happens when we read point number 2097152 - we double memory]
psxy: GMT memory: Increase 2 temporary column arrays to new length : 8388608 [happens when we read point 4194304 - we double again]
Then we finish reading and N is some number < 8388608. These arrays addresses are passed to the segment structure.
psxy: GMT memory: Initialize 2 temporary column arrays, each of length : 2097152 [OK, back to start, like at beginning of GMT we initiate new tmp arrays]
psxy: GMT memory: Increase 2 temporary column arrays to new length : 4194304 [the clipping function loop hits 2097152 and we double memory]
ERROR: Caught signal number 11 (Segmentation fault)... [but something happens before we get to the next double?

Eduardo's post: Thanks, will look.

#44 Updated by Joaquim almost 4 years ago

Recall, the deadly action occurs at

GMT_prep_tmp_arrays (GMT, j+2, 2);    /* Init or reallocate tmp vectors */
j += gmt_move_to_wesn (GMT, xx, yy, lon[i], lat[i], lon[i-1], lat[i-1], j, nx);    /* May add 2 points, which explains the j+2 stuff */

when memory is reallocated and what happens is that after the realloc xx and GMT→hidden.mem_coord0 no longer point to the same address in memory and old xx address points to an invalid memory address.

#45 Updated by Paul almost 4 years ago

Ah, I think I get it (perhaps). In gmt_wesn_clip_old we set xx, yy as shorthands to the mem_coord arrays. But I think it is true that when GMT_prep_tmp_arrays is called in that loop, it is possible that the mem_coord addresses have changed (to allocate the 128 Mb perhaps another memory chunk must be used). But xx,yy does not know that! So those point to the original address and boom.
Could you try to NOT use the xx,yy but use the full names instead, i.e

j = gmt_move_to_wesn (GMT, GMT→hidden.mem_coord[GMT_X], GMT→hidden.mem_coord[GMT_Y], lon[0], lat[0], 0.0, 0.0, 0, 0);    /* Add one point */
...
GMT→hidden.mem_coord[GMT_X][j] = xc[k];

etc, i.e., no xx, yy at all. I bet that fixes it. Sorry, off for an hour.

#46 Updated by Joaquim almost 4 years ago

Bingo. Fixed in r12859

#47 Updated by Paul almost 4 years ago

  • Status changed from Feedback to Resolved

Great, Eduardo let me know if this fixes your original issues.

#48 Updated by Joaquim almost 4 years ago

  • File eurasia.png added
  • Status changed from Resolved to In Progress

Sorry, but this not over. Right, it doesn't crash anymore but I've seen GMT produce nicer plots

#49 Updated by Paul almost 4 years ago

Perhaps possible to achieve same bad plot with a greatly reduced polygon? Then we could debug that aspect (not so easy with 4 Gb file...)

#50 Updated by Paul almost 4 years ago

  • Status changed from In Progress to Closed

Since the original bug has been fixed I am closing this issue. Please open a new issue with the clip problem on global polygons, or whatever.

Also available in: Atom PDF