Skip to content

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

@apasto

Description

@apasto

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

  • Operating system: Windows 10
  • Version of GMT: 6.4.0
  • Version of GMTMex: 2.1.0

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions