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

CPT range auto-stretching works differently with NetCDF and virtual reliefgrid inputs #5294

Closed
weiji14 opened this issue Jun 3, 2021 · 23 comments · Fixed by #5948
Closed
Assignees
Labels
bug Something isn't working

Comments

@weiji14
Copy link
Member

weiji14 commented Jun 3, 2021

Description of the problem

When comparing a grdimage plot from a regular remote file and virtual file input, while setting a specific region (e.g. Greenland/GL in the case below), the CPT range is stretched differently. Specifically, the CPT range is -3893 to 3195.5 for the remote file case, but -8592.5 to 5559 for the virtual file case. See GenericMappingTools/pygmt#750 (comment) for context.

expected (correct) actual (wrong) diff
test_grdimage_grid_and_shading_with_xarray png -expected test_grdimage_grid_and_shading_with_xarray png test_grdimage_grid_and_shading_with_xarray png -failed-diff

Need to install PyGMT branch from GenericMappingTools/pygmt#750 first using pip install https://github.com/GenericMappingTools/pygmt/archive/grdimage-xarray-shading.zip first before running the below: Edit: just do pip install pygmt=0.5.0.

Full script that generated the error

import numpy as np
import pygmt
import xarray as xr

# xarray.DataArray grid used as input for shading (-I) parameter
longitude = np.arange(0, 360, 1)
latitude = np.arange(-89, 90, 1)
x = np.sin(np.deg2rad(longitude))
y = np.linspace(start=0, stop=1, num=179)
data = y[:, np.newaxis] * x

xrgrid = xr.DataArray(
    data,
    coords=[
        ("latitude", latitude, {"units": "degrees_north"}),
        ("longitude", longitude, {"units": "degrees_east"}),
    ],
    attrs={"actual_range": [-1, 1]},
)

# Baseline (correct) case using @earth_relief_01d_g reliefgrid
fig = pygmt.Figure()
fig.grdimage(
    grid="@earth_relief_01d_g", region="GL", cmap="geo", shading=xrgrid, verbose="i"
)
fig.colorbar()
fig.show()


# Expected (wrong) case using xarray.DataArray reliefgrid
fig = pygmt.Figure()
grid = pygmt.datasets.load_earth_relief(registration="gridline")
fig.grdimage(grid=grid, region="GL", cmap="geo", shading=xrgrid, verbose="i")
fig.colorbar()
fig.show()

Full error message

Copied from https://github.com/GenericMappingTools/pygmt/runs/2733822414?check_suite_focus=true#step:11:543

Actual outcome

Shading an xarray.DataArray (via the GMT C API)

**

grdimage [INFORMATION]: Using country and state data from dcw-gmt
grdimage [INFORMATION]: Title  : DCW-GMT - The Digital Chart of the World for the Generic Mapping Tools
grdimage [INFORMATION]: Source : Processed by the GMT Team, 2018-JUL-01
grdimage [INFORMATION]: Version: 1.1.4
grdimage [INFORMATION]: Greenland
grdimage [INFORMATION]: Region implied by DCW polygons is 286.737/348.688/59.7773/83.6274
grdimage [INFORMATION]: Read intensity grid header from file @GMTAPI@-S-I-G-M-G-N-000000
grdimage [INFORMATION]: Read header from file @GMTAPI@-S-I-G-M-G-N-000001
grdimage [INFORMATION]: Spherical approximation used
grdimage [INFORMATION]: Central meridian not given, default to 317.712
grdimage [INFORMATION]: Map scale is 459.244 km per cm or 1:4.59244e+07.
grdimage [INFORMATION]: Allocate and read data from file @GMTAPI@-S-I-G-M-G-N-000001
grdimage [INFORMATION]: Importing grid data from user matrix memory location
grdimage [INFORMATION]: called with a pad < 2; skipped.
grdimage [INFORMATION]: Allocates memory and read intensity file
grdimage [INFORMATION]: Importing grid data from user matrix memory location
grdimage [INFORMATION]: Set boundary condition for all edges: natural
grdimage [INFORMATION]: Set boundary condition for left   edge: natural
grdimage [INFORMATION]: Set boundary condition for right  edge: natural
grdimage [INFORMATION]: Set boundary condition for bottom edge: natural
grdimage [INFORMATION]: Set boundary condition for top    edge: natural
grdimage [INFORMATION]: Reading CPT from File /usr/share/miniconda3/envs/pygmt/share/gmt/cpt/geo.cpt
grdimage [INFORMATION]: Auto-stretching CPT file geo to fit data range -8592.5 to 5559
grdimage [INFORMATION]: Write CPT to File /home/runner/.gmt/sessions/gmt_session.5828/gmt.78.cpt
grdimage [INFORMATION]: Save current CPT file to /home/runner/.gmt/sessions/gmt_session.5828/gmt.78.cpt !
grdimage [INFORMATION]: Evaluate image pixel colors
grdimage [INFORMATION]: Basic z(x,y) -> color image with no illumination.
grdimage [INFORMATION]: Plotting 24-bit color image
PSL: Too many colors to make colormap - using 24-bit direct color instead.
PSL: DEFLATE compressed 4725 to 4231 bytes (10.5% savings at compression level 5)

Expected outcome

Shading the @earth_relief_01d_g file directly.

grdimage [INFORMATION]: Using country and state data from dcw-gmt
grdimage [INFORMATION]: Title  : DCW-GMT - The Digital Chart of the World for the Generic Mapping Tools
grdimage [INFORMATION]: Source : Processed by the GMT Team, 2018-JUL-01
grdimage [INFORMATION]: Version: 1.1.4
grdimage [INFORMATION]: Greenland
grdimage [INFORMATION]: Region implied by DCW polygons is 286.737/348.688/59.7773/83.6274
grdimage [INFORMATION]: Read intensity grid header from file @GMTAPI@-S-I-G-M-G-N-000000
grdimage [INFORMATION]: Read header from file /home/runner/.gmt/server/earth/earth_relief/earth_relief_01d_g.grd
grdimage [INFORMATION]: Spherical approximation used
grdimage [INFORMATION]: Central meridian not given, default to 317.712
grdimage [INFORMATION]: Map scale is 459.244 km per cm or 1:4.59244e+07.
grdimage [INFORMATION]: Allocate and read data from file /home/runner/.gmt/server/earth/earth_relief/earth_relief_01d_g.grd
grdimage [INFORMATION]: Reading grid from file /home/runner/.gmt/server/earth/earth_relief/earth_relief_01d_g.grd
grdimage [INFORMATION]: All boundaries set via extended data.
grdimage [INFORMATION]: Allocates memory and read intensity file
grdimage [INFORMATION]: Importing grid data from user matrix memory location
grdimage [INFORMATION]: Set boundary condition for all edges: natural
grdimage [INFORMATION]: Set boundary condition for left   edge: natural
grdimage [INFORMATION]: Set boundary condition for right  edge: natural
grdimage [INFORMATION]: Set boundary condition for bottom edge: natural
grdimage [INFORMATION]: Set boundary condition for top    edge: natural
grdimage [INFORMATION]: Reading CPT from File /usr/share/miniconda3/envs/pygmt/share/gmt/cpt/geo.cpt
grdimage [INFORMATION]: Auto-stretching CPT file geo to fit data range -3893 to 3195.5
grdimage [INFORMATION]: Write CPT to File /home/runner/.gmt/sessions/gmt_session.5828/gmt.77.cpt
grdimage [INFORMATION]: Save current CPT file to /home/runner/.gmt/sessions/gmt_session.5828/gmt.77.cpt !
grdimage [INFORMATION]: Evaluate image pixel colors
grdimage [INFORMATION]: Basic z(x,y) -> color image with no illumination.
grdimage [INFORMATION]: Plotting 24-bit color image
PSL: Too many colors to make colormap - using 24-bit direct color instead.
PSL: DEFLATE compressed 4725 to 4375 bytes (7.4% savings at compression level 5)

Diff

-grdimage [INFORMATION]: Allocate and read data from file /home/runner/.gmt/server/earth/earth_relief/earth_relief_01d_g.grd	
-grdimage [INFORMATION]: Reading grid from file /home/runner/.gmt/server/earth/earth_relief/earth_relief_01d_g.grd	
-grdimage [INFORMATION]: All boundaries set via extended data.
+grdimage [INFORMATION]: Allocate and read data from file @GMTAPI@-S-I-G-M-G-N-000001
+grdimage [INFORMATION]: Importing grid data from user matrix memory location
+grdimage [INFORMATION]: called with a pad < 2; skipped.

-grdimage [INFORMATION]: Auto-stretching CPT file geo to fit data range -3893 to 3195.5	
-grdimage [INFORMATION]: Write CPT to File /home/runner/.gmt/sessions/gmt_session.5828/gmt.77.cpt	
-grdimage [INFORMATION]: Save current CPT file to /home/runner/.gmt/sessions/gmt_session.5828/gmt.77.cpt !
+grdimage [INFORMATION]: Auto-stretching CPT file geo to fit data range -8592.5 to 5559
+grdimage [INFORMATION]: Write CPT to File /home/runner/.gmt/sessions/gmt_session.5828/gmt.78.cpt
+grdimage [INFORMATION]: Save current CPT file to /home/runner/.gmt/sessions/gmt_session.5828/gmt.78.cpt !

System information

  • Operating system: Linux
  • GMT version (gmt --version): 6.2.0rc2
@PaulWessel
Copy link
Member

I am getting lots of these errors (the one with -1 -9 etc):

...
grdimage [ERROR]: Unrecognized option -1
grdimage [ERROR]: Unrecognized option -9
grdmix [ERROR]: Cannot find file (longitude:
grdmix [ERROR]: Cannot find file 360)>
grdimage [ERROR]: Unable to combine @GMTAPI@-S-I-G-M-G-N-000000/(longitude:/360)> into an image - aborting.
grdimage [ERROR]: Option -I: Must specify intensity file, value, or modifiers

What am I missing here?

@weiji14
Copy link
Member Author

weiji14 commented Jun 3, 2021

Ah yes, could you try doing pip install https://github.com/GenericMappingTools/pygmt/archive/grdimage-xarray-shading.zip in the terminal first. The shading feature only works on @seisman's PR branch and not on PyGMT master.

@PaulWessel
Copy link
Member

Sorry, I need more help if I am going to debug this - we are also running out on time fro 6.2.0 I think.

pip install github.com/GenericMappingTools/pygmt/archive/grdimage-xarray-shading.zip
WARNING: Requirement 'github.com/GenericMappingTools/pygmt/archive/grdimage-xarray-shading.zip' looks like a filename, but the file does not exist

I usually debug against master and now little about the pyGMT setup...

@weiji14
Copy link
Member Author

weiji14 commented Jun 3, 2021

Hmm, did you add https:// to the url?

@PaulWessel
Copy link
Member

For some reason I managed not to copy that part... OK now.

@PaulWessel
Copy link
Member

But makes no difference - same errors.

@weiji14
Copy link
Member Author

weiji14 commented Jun 3, 2021

Did you restart the Python kernel? I'm assuming you ran pip install after conda activate pygmt. Perhaps take a step back following https://github.com/GenericMappingTools/gmt/blob/master/doc/rst/source/debug.rst#debug-pygmt-in-xcode-on-macos:

conda activate pygmt  # Activate the PyGMT environment
pip install https://github.com/GenericMappingTools/pygmt/archive/grdimage-xarray-shading.zip  # install the PyGMT PR branch
python  # to get the python console
ipython  # to get interactive python console

@PaulWessel
Copy link
Member

I probably need to install ipython then

ipython
Traceback (most recent call last):
File "", line 1, in
NameError: name 'ipython' is not defined

@PaulWessel
Copy link
Member

sudo port install ipython
Password:
Error: Port ipython not found

@maxrjones
Copy link
Member

conda install ipython

@PaulWessel
Copy link
Member

OK, done, but still typing ipython into the python problem gives same error.

@weiji14
Copy link
Member Author

weiji14 commented Jun 3, 2021

Sorry, just ignore the ipython command. That's an alternative to python...

@PaulWessel
Copy link
Member

OK, getting same errors as reported at first.
is the alleged bug not reproducible with master at all?

@maxrjones
Copy link
Member

I get the bug with master gmt branch and grdimage-xarray-shading pygmt branch. I could step through to see if I notice anything if it would help.

@maxrjones
Copy link
Member

OK, getting same errors as reported at first.
is the alleged bug not reproducible with master at all?

I do not think it is reproducible with pygmt's master branch.

@weiji14
Copy link
Member Author

weiji14 commented Jun 3, 2021

OK, getting same errors as reported at first.
is the alleged bug not reproducible with master at all?

I do not think it is reproducible with pygmt's master branch.

The PR contains a feature which isn't merged into PyGMT's master branch yet, so need to install the feature branch to test.

@weiji14
Copy link
Member Author

weiji14 commented Jun 10, 2021

Ok, PR GenericMappingTools/pygmt#750 has been merged into PyGMT master so it should be easier to test things now. Note that there might be a separate bug intertwined into this, related to an offset of pixels (see GenericMappingTools/pygmt#750 (comment)).

@maxrjones
Copy link
Member

maxrjones commented Nov 3, 2021

For the example above, the gmtapi_import_grid case GMT_IS_DUPLICATE|GMT_VIA_MATRIX is used and #5940 offers a solution for the color stretching but not the lateral offset problem reported in GenericMappingTools/pygmt#750 (comment).

Here is another example that is broken in a different way:

import pygmt
# Expected (wrong) case using xarray.DataArray reliefgrid
fig = pygmt.Figure()
grid = pygmt.datasets.load_earth_relief(registration="gridline")
fig.grdimage(grid=grid, region=[0, 5, 0, 5], cmap="geo", verbose="i")
fig.colorbar()
fig.savefig("test.png")

For this example GMT_IS_REFERENCE|GMT_VIA_MATRIX is used rather than GMT_IS_DUPLICATE|GMT_VIA_MATRIX and there is no adjustment of G_obj -> header -> z_min | G_obj -> header -> z_min for the subregion. @PaulWessel, do you think this should happen in gmtlib_expand_headerpad in gmt_grdio.c, where the adjustment of the grid range and increment happens, or later when creating the colormap?

@PaulWessel
Copy link
Member

I will have a look at debugging this pygmt case.

@maxrjones
Copy link
Member

I will have a look at debugging this pygmt case.

Thanks. When I was debugging the offset issue, I used this paired down version without the shading code:

import pygmt

# Baseline (correct) case using @earth_relief_01d_g reliefgrid
fig = pygmt.Figure()
fig.grdimage(
    grid="@earth_relief_01d_g", region="GL", cmap="geo", verbose="i"
)
fig.colorbar()
fig.show()

# Expected (wrong) case using xarray.DataArray reliefgrid
fig = pygmt.Figure()
grid = pygmt.datasets.load_earth_relief(registration="gridline")
fig.grdimage(grid=grid, region="GL", cmap="geo", verbose="i")
fig.colorbar()
fig.show()

@PaulWessel
Copy link
Member

Region="GL" presumably returns the w/e/s/n in floating degrees not rounded off, yes? Could anyone try the same test with a rounded w/e/sn to the nearest integer degree and tell if the problem persists?

@weiji14
Copy link
Member Author

weiji14 commented Nov 4, 2021

Region="GL" presumably returns the w/e/s/n in floating degrees not rounded off, yes? Could anyone try the same test with a rounded w/e/sn to the nearest integer degree and tell if the problem persists?

Trying the below:

fig = pygmt.Figure()
grid = pygmt.datasets.load_earth_relief(registration="gridline")
# fig.grdimage(grid=grid, region="GL", cmap="geo", verbose="i")
# fig.grdimage(grid=grid, region="286.737/348.688/59.7773/83.6274", cmap="geo", verbose="i")
fig.grdimage(grid=grid, region="287/349/60/84", cmap="geo", verbose="i")
fig.colorbar()
fig.show()

Edit1: Sorry, I haven't tried the patch at #5940, this is just with regular GMT 6.2.

On GMT 6.2 and PyGMT v0.5.0, using integer degrees still stretches things to the global range (-8592.5 to 5559):

grdimage [INFORMATION]: Read header from file @GMTAPI@-S-I-G-M-G-N-000000
grdimage [INFORMATION]: Spherical approximation used
grdimage [INFORMATION]: Central meridian not given, default to 318
grdimage [INFORMATION]: Map scale is 459.606 km per cm or 1:4.59606e+07.
grdimage [INFORMATION]: Allocate and read data from file @GMTAPI@-S-I-G-M-G-N-000000
grdimage [INFORMATION]: Importing grid data from user matrix memory location
grdimage [INFORMATION]: called with a pad < 2; skipped.
grdimage [INFORMATION]: Reading CPT from File /srv/conda/envs/notebook/share/gmt/cpt/geo.cpt
grdimage [INFORMATION]: Auto-stretching CPT file geo to fit data range -8592.5 to 5559
grdimage [INFORMATION]: Write CPT to File /home/jovyan/.gmt/sessions/gmt_session.115/gmt.6.cpt
grdimage [INFORMATION]: Save current CPT file to /home/jovyan/.gmt/sessions/gmt_session.115/gmt.6.cpt !
grdimage [INFORMATION]: Evaluate image pixel colors
grdimage [INFORMATION]: Basic z(x,y) -> color image with no illumination.
grdimage [INFORMATION]: Plotting 24-bit color image
PSL: Too many colors to make colormap - using 24-bit direct color instead.
PSL: DEFLATE compressed 4725 to 3662 bytes (22.5% savings at compression level 5)

produces:
image

Edit2: The below are the results on GMT 6.3 from master when using integer coordinates. Range is -3824.5 to 3195.5, desired range is -3893 to 3195.5 (if using netcdf file directly):

grdimage [INFORMATION]: Read header from file @GMTAPI@-S-I-G-M-G-N-000000
grdimage [INFORMATION]: Spherical approximation used
grdimage [INFORMATION]: Central meridian not given, default to 318
grdimage [INFORMATION]: Map scale is 459.606 km per cm or 1:4.59606e+07.
grdimage [INFORMATION]: Allocate and read data from file @GMTAPI@-S-I-G-M-G-N-000000
grdimage [INFORMATION]: Subset selection for grid via a user matrix requires method GMT_IS_DUPLICATE instead of GMT_IS_REFERENCE - method has been switched
grdimage [INFORMATION]: Importing grid data from user matrix memory location
grdimage [INFORMATION]: called with a pad < 2; skipped.
grdimage [INFORMATION]: Reading CPT from File /Users/home/leongwei1/gmt-install-dir/share/cpt/geo.cpt
grdimage [INFORMATION]: Auto-stretching CPT file geo to fit data range -3824.5 to 3195.5
grdimage [INFORMATION]: Write CPT to File /Users/home/leongwei1/.gmt/sessions/gmt_session.34384/gmt.3.cpt
grdimage [INFORMATION]: Save current CPT file to /Users/home/leongwei1/.gmt/sessions/gmt_session.34384/gmt.3.cpt !
grdimage [INFORMATION]: Evaluate image pixel colors
grdimage [INFORMATION]: Basic z(x,y) -> color image with no illumination.
grdimage [INFORMATION]: Plotting 24-bit color image
PSL: Too many colors to make colormap - using 24-bit direct color instead.
PSL: DEFLATE compressed 4725 to 3737 bytes (20.9% savings at compression level 5)

produces

image

@maxrjones
Copy link
Member

In the example below, the colormap is not correct for virtual files even after the last two PRs:

import pygmt

# Baseline (correct) case using @earth_relief_01d_g reliefgrid
fig = pygmt.Figure()
fig.grdimage(
    grid="@earth_relief_01d_g", region=[0, 5, 0, 5], cmap="geo", verbose="i"
)
fig.colorbar()
fig.savefig("test_correct.png")

# Expected (wrong) case using xarray.DataArray reliefgrid
fig = pygmt.Figure()
grid = pygmt.datasets.load_earth_relief(registration="gridline")
fig.grdimage(grid=grid, region=[0, 5, 0, 5], cmap="geo", verbose="i")
fig.colorbar()
fig.savefig("test_wrong.png")

test_correct.png:

test_wrong.png:

In constrast to the original example, the incorrect case here goes through the GMT_IS_REFERENCE|GMT_VIA_MATRIX case:

gmt/src/gmt_api.c

Lines 5242 to 5327 in 8c55a80

case GMT_IS_REFERENCE|GMT_VIA_MATRIX: /* The user's 2-D grid array of some sort, + info in the args [NOT YET FULLY TESTED] */
/* Getting a matrix info S_obj->resource. Create grid header and then pass the grid pointer via the matrix pointer */
if ((M_obj = S_obj->resource) == NULL) return_null (API, GMT_PTR_IS_NULL);
/* If passed by reference and the module requires a pad then we must switch to duplication */
if (mode & GMT_GRID_NEEDS_PAD1 || mode & GMT_GRID_NEEDS_PAD2) { /* Cannot do this by reference, switch to duplication */
start_over_method = gmtapi_switch_method (API, S_obj, &method, "Grid via a user matrix requires method GMT_IS_DUPLICATE instead of GMT_IS_REFERENCE due to a padding requirement - method has been switched\n");
goto start_over_import_grid;
}
/* This method requires the input data to be a GMT_GRD_FORMAT matrix - otherwise we should be DUPLICATING */
if (!(M_obj->shape == GMT_IS_ROW_FORMAT && M_obj->type == GMT_GRDFLOAT && (mode & GMT_GRID_IS_COMPLEX_MASK) == 0)) {
start_over_method = gmtapi_switch_method (API, S_obj, &method, "Grid via a user matrix requires method GMT_IS_DUPLICATE instead of GMT_IS_REFERENCE due to incompatible data type for a grid - method has been switched\n");
goto start_over_import_grid;
}
/* Determine if it is possible to use the matrix given the region selected and the fact we chose GMT_IS_REFERENCE. This test will
* only kick in after we allocate the G_obj and come back the second time (after getting header) since otherwise S_obj->wesn is not set yet */
if (!(!S_obj->region ||
(S_obj->wesn[XLO] >= M_obj->range[XLO] && S_obj->wesn[XHI] <= M_obj->range[XHI] && S_obj->wesn[YLO] >= M_obj->range[YLO] && S_obj->wesn[YHI] <= M_obj->range[YHI]) ||
gmt_whole_earth (GMT, M_obj->range, S_obj->wesn))) { /* Cannot do this by reference, switch to duplication */
start_over_method = gmtapi_switch_method (API, S_obj, &method, "Subset selection for grid via a user matrix requires method GMT_IS_DUPLICATE instead of GMT_IS_REFERENCE - method has been switched\n");
goto start_over_import_grid;
}
if (grid == NULL) { /* Only allocate when not already allocated, and only get container. Note: Cannot have pad since input matrix won't have one */
uint64_t dim[3] = {M_obj->n_rows, M_obj->n_columns, 1};
if ((G_obj = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, dim, M_obj->range, M_obj->inc, M_obj->registration, 0, NULL)) == NULL)
return_null (API, GMT_MEMORY_ERROR);
}
else
G_obj = grid;
HH = gmt_get_H_hidden (G_obj->header);
G_obj->header->complex_mode = (mode & GMT_GRID_IS_COMPLEX_MASK); /* Set the complex mode */
done = (mode & GMT_CONTAINER_ONLY) ? false : true; /* Not done until we read grid */
if (! (mode & GMT_DATA_ONLY)) {
gmtapi_matrixinfo_to_grdheader (GMT, G_obj->header, M_obj); /* Populate a GRD header structure */
if (HH->grdtype > GMT_GRID_CARTESIAN) gmt_set_geographic (GMT, GMT_IN);
/* Temporarily set data pointer for convenience; removed later */
#ifdef DOUBLE_PRECISION_GRID
G_obj->data = M_obj->data.f8;
#else
G_obj->data = M_obj->data.f4;
#endif
G_obj->header->z_min = +DBL_MAX;
G_obj->header->z_max = -DBL_MAX;
HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
gmt_M_grd_loop (GMT, G_obj, row, col, ij) {
if (gmt_M_is_fnan (G_obj->data[ij]))
HH->has_NaNs = GMT_GRID_HAS_NANS;
else {
G_obj->header->z_min = MIN (G_obj->header->z_min, G_obj->data[ij]);
G_obj->header->z_max = MAX (G_obj->header->z_max, G_obj->data[ij]);
}
}
G_obj->data = NULL; /* Since data are not requested yet */
if (mode & GMT_CONTAINER_ONLY) /* Just needed the header but had to set zmin/zmax first */
break;
}
if ((new_ID = gmtapi_get_object (API, GMT_IS_GRID, G_obj)) == GMT_NOTSET)
return_null (API, GMT_OBJECT_NOT_FOUND);
if ((new_item = gmtlib_validate_id (API, GMT_IS_GRID, new_ID, GMT_IN, GMT_NOTSET)) == GMT_NOTSET)
return_null (API, GMT_OBJECT_NOT_FOUND);
GMT_Report (API, GMT_MSG_INFORMATION, "Referencing grid data from user memory location\n");
#ifdef DOUBLE_PRECISION_GRID
G_obj->data = M_obj->data.f8;
#else
G_obj->data = M_obj->data.f4;
#endif
GH = gmt_get_G_hidden (G_obj);
MH = gmt_get_M_hidden (M_obj);
S_obj->alloc_mode = MH->alloc_mode; /* Pass on alloc_mode of matrix */
GH->alloc_mode = GMT_ALLOC_EXTERNALLY; /* Since we cannot have both M and G try to free */
API->object[new_item]->resource = G_obj;
API->object[new_item]->status = GMT_IS_USED; /* Mark as read */
GH->alloc_level = API->object[new_item]->alloc_level; /* Since allocated here */
if (gmt_whole_earth (GMT, M_obj->range, S_obj->wesn)) {
/* Global grids passed via matrix are not rotated to fit the desired global region, so we need to correct the wesn for this grid to match the matrix */
gmt_M_memcpy (G_obj->header->wesn, M_obj->range, 4U, double);
}
else if (S_obj->region) { /* Possibly adjust the pad so inner region matches wesn */
if (S_obj->reset_pad) { /* First undo a prior sub-region used with this memory grid */
gmtlib_contract_headerpad (GMT, G_obj->header, S_obj->orig_pad, S_obj->orig_wesn);
S_obj->reset_pad = 0;
}
if (gmtlib_expand_headerpad (GMT, G_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn))
S_obj->reset_pad = 1;
}
break;

The other header attributes are updated in gmtlib_expand_headerpad, so I think that the header z_min and z_max could get updated here as well.

gmt/src/gmt_grdio.c

Lines 3235 to 3264 in 8c55a80

unsigned int gmtlib_expand_headerpad (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, double *new_wesn, unsigned int *orig_pad, double *orig_wesn) {
unsigned int tmp_pad[4] = {0, 0, 0, 0}, delta[4] = {0, 0, 0, 0}, k = 0;
/* When using subset with memory grids we cannot actually cut the grid but instead
* must temporarily change the pad to match the desired inner region wesn. This means
* the pads will change and can be quite large. */
struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h);
gmt_M_memcpy (tmp_pad, h->pad, 4, unsigned int); /* Initialize new pad to the original pad */
/* First determine which (and how many, k) of the 4 new boundaries are inside the original region and update the padding: */
if (new_wesn[XLO] > h->wesn[XLO]) k++, tmp_pad[XLO] += urint ((new_wesn[XLO] - h->wesn[XLO]) * HH->r_inc[GMT_X]);
if (new_wesn[XHI] < h->wesn[XHI]) k++, tmp_pad[XHI] += urint ((h->wesn[XHI] - new_wesn[XHI]) * HH->r_inc[GMT_X]);
if (new_wesn[YLO] > h->wesn[YLO]) k++, tmp_pad[YLO] += urint ((new_wesn[YLO] - h->wesn[YLO]) * HH->r_inc[GMT_Y]);
if (new_wesn[YHI] < h->wesn[YHI]) k++, tmp_pad[YHI] += urint ((h->wesn[YHI] - new_wesn[YHI]) * HH->r_inc[GMT_Y]);
if (k) { /* Yes, pad will change since region is different for k of the 4 sides */
for (k = 0; k < 4; k++) delta[k] = tmp_pad[k] - h->pad[k]; /* Columns with data being passed as padding */
gmt_M_memcpy (orig_pad, h->pad, 4, unsigned int); /* Place the original grid pad in the provided array */
gmt_M_memcpy (orig_wesn, h->wesn, 4, double); /* Place the original grid wesn in the provided array */
gmt_M_memcpy (h->pad, tmp_pad, 4, unsigned int); /* Place the new pad in the grid header */
gmt_M_memcpy (h->wesn, new_wesn, 4, double); /* Place the new wesn in the grid header */
gmt_set_grddim (GMT, h); /* This recomputes n_columns|n_rows. */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtlib_expand_headerpad: %d pad sides changed. Now %u/%u/%u/%u\n",
k, h->pad[XLO], h->pad[XHI], h->pad[YLO], h->pad[YHI]);
for (k = 0; k < 4; k++) { /* If pad now contains data then change the BC to reflect this */
if (delta[k] >= orig_pad[k]) HH->BC[k] = GMT_BC_IS_DATA;
}
}
else
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtlib_expand_headerpad: No pad adjustment needed\n");
return k;
}

PaulWessel added a commit that referenced this issue Nov 4, 2021
When a matrix is passed by reference and there is no need to switch to duplication, we adjust the grid's pad to simulate a cut.  However, we never updated the header zmin/zmax to reflect the new data region.  This PR takes care of this for grids (and cubes).  Closes #5294.
PaulWessel added a commit that referenced this issue Nov 4, 2021
* Update the z min/max inside readonly subsets

When a matrix is passed by reference and there is no need to switch to duplication, we adjust the grid's pad to simulate a cut.  However, we never updated the header zmin/zmax to reflect the new data region.  This PR takes care of this for grids (and cubes).  Closes #5294.

* use GMT_LOCAL
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.

3 participants