search

GenericMappingTools / gmt

The Generic Mapping Tools
https://www.generic-mapping-tools.org
Other
843 stars 352 forks source link

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

Closed seisman closed 2 months ago

seisman commented 3 months ago

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 commented 3 months ago

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 commented 3 months ago

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 commented 3 months ago

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 commented 3 months ago

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 commented 3 months ago

A yes, the pad.

joa-quim commented 3 months ago

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 commented 3 months ago

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 commented 2 months ago

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 commented 2 months ago

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

seisman commented 2 months ago

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

seisman commented 2 months ago

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: https://github.com/GenericMappingTools/gmt/blob/206db56471435d75825dd72869d63b26664fe20a/src/grdinfo.c#L823-L840

The min/max in the grid header is only updated when -M option is set: https://github.com/GenericMappingTools/gmt/blob/206db56471435d75825dd72869d63b26664fe20a/src/grdinfo.c#L855

When output the grid information, grdinfo always uses the information from the header: https://github.com/GenericMappingTools/gmt/blob/206db56471435d75825dd72869d63b26664fe20a/src/grdinfo.c#L967

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 commented 2 months ago

I agrre. Line 855 should become

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