grdinfo -L0 doesn't work for external wrappers about gmt

[joa-quim]

Hmm, it works for me in Julia

julia> G = mat2grid(reshape(collect(0.0f0:8*9-1), 8,9))
A GMTgrid object with 1 layers of type Float32
Gridline node registration used
x_min: 1.0      x_max :9.0      x_inc :1.0      n_columns :9
y_min: 1.0      y_max :8.0      y_inc :1.0      n_rows :8
z_min: 0.0      z_max :71.0
Mem layout:     BCB

julia> grdinfo(G, L=0, C=true)

1×12 GMTdataset{Float64, 2}
 Row │ x_min  x_max  y_min  y_max  z_min  z_max   dx   dy  n_cols  n_rows  reg  isgeog
─────┼─────────────────────────────────────────────────────────────────────────────────
   1 │   1.0    9.0    1.0    8.0    0.0   71.0  1.0  1.0     9.0     8.0  0.0     0.0


julia> grdinfo(G,  C=true)

1×12 GMTdataset{Float64, 2}
 Row │ x_min  x_max  y_min  y_max  z_min  z_max   dx   dy  n_cols  n_rows  reg  isgeog
─────┼─────────────────────────────────────────────────────────────────────────────────
   1 │   1.0    9.0    1.0    8.0    0.0   71.0  1.0  1.0     9.0     8.0  0.0     0.0

uint64_t dim[2] = {5, 4};

Shouldn't this be {8, 9} (or {9, 8}), the size of the grid?

[seisman]

Shouldn't this be {8, 9} (or {9, 8}), the size of the grid?

You're right. It's a typo in the original test.

After changing to dim[2] = {9, 8}, the result is still incorrect:

0	0	0	0	0	0	0	0	1	1	0	0

[seisman]

If I change GMT_IN to GMT_IN|GMT_IS_REFERENCE, the result is slightly better:

0	8	0	7	NaN	NaN	1	1	9	8	0	0

[seisman]

Shouldn't this be {8, 9} (or {9, 8}), the size of the grid?

You're right. It's a typo in the original test.

After changing to dim[2] = {9, 8}, the result is still incorrect:

0	0	0	0	0	0	0	0	1	1	0	0

Actually, {5, 4} might be the correct size without counting the padding.

For example, here is the original testapi_gmtgrid.c file in the src directory:

#include "gmt.h"
/*
 * Testing the use of user grid memory allocated externally.
 */

int main () {
	unsigned int mode = GMT_SESSION_EXTERNAL;
	struct GMT_GRID *G = NULL;
	struct GMTAPI_CTRL *API = NULL;
	uint64_t dim[2] = {5, 4};
	char args[1000] = {""};
	char input[GMT_VF_LEN] = {""};
	gmt_grdfloat *data = NULL;

	API = GMT_Create_Session ("testapi_gmtgrid", 2U, mode, NULL);

	/* Create a container for the grid */
	G = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, dim, NULL, NULL, 0, 2, NULL);

	/* Allocate memory of the data array in the external program (C or PyGMT) */
	data = (gmt_grdfloat *)malloc(sizeof(gmt_grdfloat) * 9 * 8);
	for (int i = 0; i < 9 * 8; i++) data[i] = (gmt_grdfloat)i;

	/* Assign the user data to the GMT_GRID structure */
	G->data = data;

	/* Using the grid */
	GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN, G, input);
	sprintf (args, "%s", input);
	GMT_Call_Module (API, "grd2xyz", GMT_MODULE_CMD, args);
	GMT_Close_VirtualFile (API, input);

	/* Free the data array allocated by the external program */
	free (data);

	if (GMT_Destroy_Session (API)) return EXIT_FAILURE;
	exit (0);
}

Its output is:

0	3	20
1	3	21
2	3	22
3	3	23
4	3	24
0	2	29
1	2	30
2	2	31
3	2	32
4	2	33
0	1	38
1	1	39
2	1	40
3	1	41
4	1	42
0	0	47
1	0	48
2	0	49
3	0	50
4	0	51

After changing dim[2] to {9, 8}, the output is:

0	7	28
1	7	29
2	7	30
3	7	31
4	7	32
5	7	33
6	7	34
7	7	35
8	7	36
0	6	41
1	6	42
2	6	43
3	6	44
4	6	45
5	6	46
6	6	47
7	6	48
8	6	49
0	5	54
1	5	55
2	5	56
3	5	57
4	5	58
5	5	59
6	5	60
7	5	61
8	5	62
0	4	67
1	4	68
2	4	69
3	4	70
4	4	71
5	4	-1.47249807766e+19
6	4	4.20389539297e-45
7	4	1.1350517561e-43
8	4	0
0	3	0
1	3	0
2	3	4.63080422399e+27
3	3	7.17755831329e+22
4	3	2.73762157729e+20
5	3	2.96097025285e+29
6	3	208860992
7	3	2.73446536319e+20
8	3	4.54482212584e+30
0	2	0
1	2	7.41286887628e-43
2	2	0
3	2	1.25415113939e-37
4	2	0
5	2	-2.47539658871e+35
6	2	1.78646230296e-38
7	2	0
8	2	0
0	1	0
1	1	0
2	1	0
3	1	0
4	1	0
5	1	0
6	1	0
7	1	0
8	1	0
0	0	0
1	0	0
2	0	0
3	0	0
4	0	0
5	0	0
6	0	0
7	0	0
8	0	0

[joa-quim]

A yes, the pad.

[joa-quim]

Once more: yes, the pad

The example is correct without the -L0 because it only prints what's in the grid's header. But with it, it scans the grid and that's where the problem arises. With a pad = 0 it actually works well too, but with pad = 2 in gmt_api.c#L5383

if (!S_obj->region && gmt_grd_pad_status (GMT, G_obj->header, GMT->current.io.pad)) {	/* Want an exact copy with no subset and same padding */

GMT->current.io.pad is {0,0,0,0} when it should be {2,2,2,2} and the code flow, flows along a deadly branch, the wesn is set to [0,0,0,0] and grid (from virtual file) is wrongly read.

Now, is this a flaw of this simple example or is it a more fundamental issue in the API? I almost don't use the xx_VirtualFile in Julia so cannot tell much about it.

[joa-quim]

Or maybe the problem is in grdinfo line 692 that always sets the session pad to 0

gmt_set_pad (API->GMT, 0);	/* Not read pads to simplify operations */

[seisman]

Or maybe the problem is in grdinfo line 692 that always sets the session pad to 0

gmt_set_pad (API->GMT, 0);	/* Not read pads to simplify operations */

I tried to comment out this line, and the new result is:

0	4	0	3	NaN	NaN	1	1	5	4	0	0

[joa-quim]

Which is the correct result. The NaNs originate in the test example that does not set the grid's z_min|max

[seisman]

But -L0 means Report range of data by actually reading them (not from header)., So min/max should be correctly reported, no?

[seisman]

Looking at the grdinfo source code.

grdinfo scans all the grid nodes to find the min/max values of the grid and the min/max values are stored in the variables v_min/v_max:

gmt/src/grdinfo.c

Lines 823 to 840 in 206db56

	v_min = DBL_MAX; v_max = -DBL_MAX;
	ij_min = ij_max = n = 0;
	for (level = 0; level < header->n_bands; level++) {
	for (row = 0; row < header->n_rows; row++) {
	for (col = 0, ij = gmt_M_ijp (header, row, 0); col < header->n_columns; col++, ij++) {
	node = ij + here;
	if (gmt_M_is_fnan (data[node])) continue;
	if (data[node] < v_min) {
	v_min = data[node]; ij_min = ij; if (is_cube) z_min = U->z[level];
	}
	if (data[node] > v_max) {
	v_max = data[node]; ij_max = ij; if (is_cube) z_max = U->z[level];
	}
	n++;
	}
	}
	here += header->size;
	}

The min/max in the grid header is only updated when -M option is set:

gmt/src/grdinfo.c

Line 855 in 206db56

if (Ctrl->M.mode == GRDINFO_FORCE) { /* Update header */

When output the grid information, grdinfo always uses the information from the header:

gmt/src/grdinfo.c

Line 967 in 206db56

out[col++] = header->z_min; out[col++] = header->z_max;

I think this is a bug. When -L0 is given, we should do either of these:

1. Update the header->z_min/header->z_max with v_min/v_max
2. Output v_min/v_max rather than header->z_min/header->z_max