search

GenericMappingTools / gmtmex

GMT API for MATLAB
https://github.com/GenericMappingTools/gmtmex/wiki
GNU Lesser General Public License v2.1
42 stars 18 forks source link

grdcut -N on larger region results in no grid extension #31

Open apasto opened 2 years ago

apasto commented 2 years ago

Hello everyone,

Description of the problem

I am trying to extend a grid onto a new region, including but exceeding the original boundaries, using grdcut with the -N option. However, the returned grid structure retains the original region unchanged: no nan-padding occurs.

Full script that generated the error

This minimal example reproduces the issue:

region_0 = [-5, 5, -10, 10];
n_el = [21, 41];

x = linspace(region_0(1), region_0(2), n_el(1));
y = linspace(region_0(3), region_0(4), n_el(2));
inc = [x(2) - x(1), y(2) - y(1)];
z = ones(length(y), length(x));
header = [region_0, min(z(:)), max(z(:)), 0, inc];

G_0 = gmt('wrapgrid', z, header);

region_1 = [-6, 6, -11, 11];
G_1 = gmt(['grdcut -Vd -N -R', sprintf('%g/%g/%g/%g', region_1)], G_0);

assert(all(region_1 == G_1.range(1:4)))

This also fails with fill values other than nan, e.g. -N0.

G_1 region is left unchanged, still equal to G_0 extents:

>> G_0.range
ans =    -5     5   -10    10     1     1
>> G_1.range
ans =    -5     5   -10    10     1     1

Writing G_0 to file and calling gmt grdcut -N outside Matlab, from the shell, the nan-fill extension works as expected.

Full error message

grdcut [DEBUG]: Got regular w/e/s/n for region (-6/6/-11/11)
grdcut [INFORMATION]: Processing input grid
grdcut [DEBUG]: Set_Object for family: 1
grdcut [INFORMATION]: Cartesian input grid
grdcut [INFORMATION]: Requested subset exceeds data domain on the left side - nodes in the extra area will be initialized to -nan(ind)
grdcut [INFORMATION]: Requested subset exceeds data domain on the right side - nodes in the extra area will be initialized to -nan(ind)
grdcut [INFORMATION]: Requested subset exceeds data domain on the bottom side - nodes in the extra area will be initialized to -nan(ind)
grdcut [INFORMATION]: Requested subset exceeds data domain on the top side - nodes in the extra area will be initialized to -nan(ind)
grdcut [DEBUG]: Set_Object for family: 1
grdcut [INFORMATION]: Cartesian input grid
grdcut [DEBUG]: gmtapi_expand_headerpad: No pad adjustment needed
grdcut [DEBUG]: Chosen boundary condition for all edges: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for all edges: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for left   edge: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for right  edge: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for bottom edge: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for top    edge: natural
grdcut [DEBUG]: Object ID 17 : Registered Grid Memory Reference f3c22b56c0 as an Input resource with geometry Surface [n_objects = 3]
grdcut [DEBUG]: Successfully duplicated a Grid
==> 3 API Objects at end of GMT_Duplicate_Data
--------------------------------------------------------
K.. ID RESOURCE.... FAMILY.... ACTUAL.... DIR... S O M L
--------------------------------------------------------
* 0 15   f3c22bfe30 Grid       Grid       Input  0 Y N 0
* 1 16   f3c22bfbc0 Grid       Grid       Output 0 Y Y 0
* 2 17   f3c22b56c0 Grid       Grid       Input  0 Y N 1
--------------------------------------------------------
grdcut [INFORMATION]: File spec:    W E S N dx dy n_columns n_rows:
grdcut [INFORMATION]: Old:grdcut [INFORMATION]:     -5  5   -10 10  0.5 0.5 21  41
grdcut [INFORMATION]: New:grdcut [INFORMATION]:     -6  6   -11 11  0.5 0.5 25  45
grdcut [DEBUG]: GMT_Write_Data: Writing Grid to memory object 16 from object 17 which transfers ownership
grdcut [DEBUG]: gmtapi_begin_io: Output resource access is now enabled [container]
grdcut [DEBUG]: gmtapi_export_data: Messenger dummy output container for object 16 [item 1] freed and set resource=data=NULL
grdcut [DEBUG]: gmtapi_export_grid: Passed ID = 16 and mode = 0
grdcut [INFORMATION]: Referencing grid data to GMT_GRID memory location
grdcut [DEBUG]: Chosen boundary condition for all edges: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for all edges: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for left   edge: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for right  edge: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for bottom edge: natural
grdcut [INFORMATION]: gmt_grd_BC_set: Set boundary condition for top    edge: natural
grdcut [DEBUG]: GMT_End_IO: Output resource access is now disabled
grdcut [DEBUG]: Set_Object for family: 1
==> 3 API Objects at end of GMT_Write_Data
--------------------------------------------------------
K.. ID RESOURCE.... FAMILY.... ACTUAL.... DIR... S O M L
--------------------------------------------------------
* 0 15   f3c22bfe30 Grid       Grid       Input  0 Y N 0
* 1 16   f3c22b56c0 Grid       Grid       Output 2 Y N 0
* 2 17   f3c22b56c0 Grid       Grid       Input  0 N N 1
--------------------------------------------------------
==> 3 API Objects at end of GMTAPI_Garbage_Collection entry
--------------------------------------------------------
K.. ID RESOURCE.... FAMILY.... ACTUAL.... DIR... S O M L
--------------------------------------------------------
* 0 15   f3c22bfe30 Grid       Grid       Input  0 Y N 0
* 1 16   f3c22b56c0 Grid       Grid       Output 2 Y N 0
* 2 17   f3c22b56c0 Grid       Grid       Input  0 N N 1
--------------------------------------------------------
grdcut [DEBUG]: gmtlib_unregister_io: Unregistering object no 17 [n_objects = 2]
==> 2 API Objects at end of GMTAPI_Garbage_Collection exit
--------------------------------------------------------
K.. ID RESOURCE.... FAMILY.... ACTUAL.... DIR... S O M L
--------------------------------------------------------
* 0 15   f3c22bfe30 Grid       Grid       Input  0 Y N 0
* 1 16   f3c22b56c0 Grid       Grid       Output 2 Y N 0
--------------------------------------------------------

System information

welcome[bot] commented 2 years ago

👋 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!

joa-quim commented 2 years ago

Yep, can confirm this behavior in Julia (GMT.jl) too, which functions in a similar way to the gmtmex. This points to the issue being in the GMT lib itself.

Thanks for the report.