Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

grdinfo -L doesn't work for external wrappers #8525

Closed
seisman opened this issue Jun 20, 2024 · 12 comments · Fixed by #8528
Closed

grdinfo -L doesn't work for external wrappers #8525

seisman opened this issue Jun 20, 2024 · 12 comments · Fixed by #8528
Labels
bug Something isn't working
Milestone

Comments

@seisman
Copy link
Member

seisman commented Jun 20, 2024

Here is the C source codes that is modified from the testapi_gmtgrid.c function:

#include "gmt.h"

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 -L0 -C", input);
	GMT_Call_Module (API, "grdinfo", 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);
}

Compiling and running the C program returns:

0	0	0	0	47	47	0	0	1	1	0	0

Without -L0, it returns:

0	4	0	3	NaN	NaN	1	1	5	4	0	0
@joa-quim
Copy link
Member

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
Copy link
Member Author

seisman commented Jun 21, 2024

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
Copy link
Member Author

seisman commented Jun 21, 2024

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
Copy link
Member Author

seisman commented Jun 21, 2024

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
Copy link
Member

A yes, the pad.

@joa-quim
Copy link
Member

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
Copy link
Member

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
Copy link
Member Author

seisman commented Jun 24, 2024

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
Copy link
Member

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

@seisman
Copy link
Member Author

seisman commented Jun 24, 2024

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

@seisman
Copy link
Member Author

seisman commented Jun 24, 2024

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:

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

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

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

@joa-quim
Copy link
Member

I agrre. Line 855 should become

if (Ctrl->M.mode == GRDINFO_FORCE || Ctrl->L.active)

@seisman seisman added this to the 6.6.0 milestone Sep 30, 2024
@seisman seisman changed the title grdinfo -L0 doesn't work for external wrappers grdinfo -L doesn't work for external wrappers Dec 4, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants