grdfill fails, probably in situation where holes extend to the edge of the grid domain

[mlees0209]

I have found that grdfill is failing, running into a segmentation fault. I am running gmt 6.4.0 on MacOS. It fails when trying to do certain operations on a grid which has lots of nans running up against the edge of the grid. It fails when using -L and -Ac but not -An. (I have not tried other options)

Command that creates failure: (the file D20191226.nc can be found here: https://drive.google.com/file/d/1rDbvuEidaeIdsuik0O7kSReNcsjzoVt2/view?usp=sharing)

gmt grdfill D20191226.nc -Ac -GD20191226_HOLESFILLED.nc -Vd

Full error message

gmt [DEBUG]: GMT_Create_Session: Terminal width = 122 gmt [DEBUG]: Obtained the ppid from parent: 36863 gmt [DEBUG]: Enter: gmtinit_new_GMT_ctrl gmt [DEBUG]: GMT->session.SHAREDIR = /Users/mlees/.conda/envs/general_science/share/gmt gmt [DEBUG]: GMT->session.HOMEDIR = /Users/mlees gmt [DEBUG]: GMT->session.USERDIR = /Users/mlees/.gmt [created] gmt [DEBUG]: GMT->session.CACHEDIR = /Users/mlees/.gmt/cache [created] gmt [DEBUG]: GMT: 0. Will try to find subdir=postscriptlight stem = PSL_custom_fonts suffix=.txt gmt [DEBUG]: GMT: 1. gmt_getsharepath trying current dir gmt [DEBUG]: GMT: 2. gmt_getsharepath trying USERDIR /Users/mlees/.gmt gmt [DEBUG]: GMT: 3. gmt_getsharepath trying USERDIR subdir /Users/mlees/.gmt/postscriptlight gmt [DEBUG]: GMT: 4. gmt_getsharepath trying SHAREDIR subdir /Users/mlees/.conda/envs/general_science/share/gmt/postscriptlight gmt [DEBUG]: GMT: 5. gmt_getsharepath trying SHAREDIR /Users/mlees/.conda/envs/general_science/share/gmt gmt [DEBUG]: GMT: 6. gmt_getsharepath failed gmt [DEBUG]: Map distance calculation will be Cartesian gmt [DEBUG]: Exit: gmtinit_new_GMT_ctrl gmt [DEBUG]: Enter: New_PSL_Ctrl gmt [DEBUG]: Exit: New_PSL_Ctrl gmt [DEBUG]: Enter: gmt_manage_workflow gmt [DEBUG]: Exit : gmt_manage_workflow gmt [DEBUG]: Enter: PSL_beginsession gmt [DEBUG]: Exit : PSL_beginsession gmt [DEBUG]: Enter: PSL_setdefaults gmt [DEBUG]: Exit : PSL_setdefaults gmt [DEBUG]: Enter: gmtlib_io_init gmt [DEBUG]: Exit : gmtlib_io_init gmt [DEBUG]: Enter: gmt_hash_init gmt [DEBUG]: Exit: gmt_hash_init gmt [DEBUG]: Enter: gmt_hash_init gmt [DEBUG]: Exit: gmt_hash_init gmt [DEBUG]: Enter: gmt_reload_settings gmt [DEBUG]: The PROJ_GEODESIC set to Vincenty gmt [DEBUG]: Look for file /Users/mlees/gmt.conf gmt [DEBUG]: Look for file /Users/mlees/.gmt/gmt.conf gmt [DEBUG]: Look for file /Users/mlees/.gmt/server/gmt.conf gmt [DEBUG]: Look for file /Users/mlees/.gmt/cache/gmt.conf gmt [DEBUG]: Could not find file gmt.conf gmt [DEBUG]: Exit: gmt_reload_settings gmt [DEBUG]: Enter: gmtlib_plot_C_format gmt [DEBUG]: Exit: gmtlib_plot_C_format gmt [DEBUG]: Enter: gmtinit_get_history gmt [DEBUG]: Enter: gmt_hash_init gmt [DEBUG]: Exit: gmt_hash_init gmt [DEBUG]: Exit: gmtinit_get_history gmt [DEBUG]: Initialize FFTW with 8 threads. gmt [DEBUG]: GMT_Create_Session initialized GMT structure gmt [DEBUG]: Loading core GMT shared library: libgmt.dylib gmt [DEBUG]: Shared Library # 0 (core). Path = libgmt.dylib gmt [DEBUG]: Loading GMT plugins from: /Users/mlees/.conda/envs/general_science/lib/gmt/plugins gmt [DEBUG]: Shared Library # 1 (supplements). Path = /Users/mlees/.conda/envs/general_science/lib/gmt/plugins/supplements.so gmt [DEBUG]: Revised options: D20191226.nc -Ac0.0 -GD20191226_HOLESFILLED.nc -Vd grdfill [DEBUG]: Found readable file D20191226.nc grdfill [DEBUG]: Replace file D20191226.nc with path D20191226.nc grdfill [DEBUG]: Replace file D20191226.nc with D20191226.nc grdfill [DEBUG]: Found readable file D20191226.nc grdfill [DEBUG]: Replace file D20191226.nc with path D20191226.nc grdfill [DEBUG]: Found readable file D20191226.nc grdfill [DEBUG]: Replace file D20191226.nc with path D20191226.nc grdfill [DEBUG]: Found readable file D20191226.nc grdfill [DEBUG]: Object ID 0 : Registered Grid File D20191226.nc as an Input resource with geometry Surface [n_objects = 1] grdfill [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container] grdfill [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 1 grdfill [DEBUG]: Found readable file D20191226.nc grdfill [DEBUG]: Calling nc_open on D20191226.nc, ncid = 65536, err = 0 grdfill [DEBUG]: Calling nc_close on ncid 65536, err = 0 grdfill [DEBUG]: Calling nc_open on D20191226.nc, ncid = 65536, err = 0 grdfill [DEBUG]: Calling nc_close on ncid 65536, err = 0 grdfill [DEBUG]: Call gmtgrdio_doctor_geo_increments on a geographic grid grdfill [DEBUG]: Geographic input grid, longitudes span less than 360 grdfill [DEBUG]: GMT_End_IO: Input resource access is now disabled grdfill [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container] grdfill [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 2 grdfill [INFORMATION]: Reading grid from file D20191226.nc grdfill [DEBUG]: Calling nc_open on D20191226.nc, ncid = 65536, err = 0 grdfill [DEBUG]: packed z-range: [-0.179119,0.0716667] grdfill [DEBUG]: Calling nc_close on ncid 65536, err = 0 grdfill [DEBUG]: Geographic input grid, longitudes span less than 360 grdfill [DEBUG]: Chosen boundary condition for all edges: geographic grdfill [INFORMATION]: gmt_grd_BC_set: Set boundary condition for all edges: natural grdfill [INFORMATION]: gmt_grd_BC_set: Set boundary condition for left edge: natural grdfill [INFORMATION]: gmt_grd_BC_set: Set boundary condition for right edge: natural grdfill [INFORMATION]: gmt_grd_BC_set: Set boundary condition for bottom edge: natural grdfill [INFORMATION]: gmt_grd_BC_set: Set boundary condition for top edge: natural grdfill [DEBUG]: GMT_End_IO: Input resource access is now disabled Segmentation fault: 11

[welcome[bot]]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

[PaulWessel]

Not a bug per se, but your grid has tons of NaNs all around a smaller feature, and what happens is that we use a recursive John Von Neumann neighbourhood algorithm and it ends up calling itself tons of time until some recursive stack is exhausted.

By resampling your grid down for testing it worked fine. Other than us rewriting the algorithm without recursing my workaround solutions would be

1. Use grdcut to extract the subset that skips all the boundary NaNs
2. Run grdfill on this smaller grid
3. Extend the filled grid with what ever you think is sensible in the large area theatre all NaN.

[mlees0209]

Thanks for the speedy response! I found another workaround which was to use grdclip to replace all the NaNs with zeros instead.

I suppose I was using the wrong tool for the matter at hand. I did not know that grdcut can extract a subset which removes all the boundary NaNs. That's helpful.

Could an improved error message or documentation be developed to guide others away from grdfill when the grid has a very large number of NaNs?

[PaulWessel]

Will see if possible to know the recursive limit

[joa-quim]

Can the Tail-Call optimization technique explained in this ChatGPT post(?!) be used to make the stack size problem go away? Basically what it says is that the grdfill_trace_the_hole should be re-written to make the last line a call to itself and the test to stop the recursion be in the function body.

Also of interest this SO thread

[joa-quim]

Another occurrence https://forum.generic-mapping-tools.org/t/segmentation-fault-in-grdfill/5037