search

GenericMappingTools / gmt

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

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

Closed weiji14 closed 2 years ago

weiji14 commented 3 years ago

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 https://github.com/GenericMappingTools/pygmt/pull/750#discussion_r644475505 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 https://github.com/GenericMappingTools/pygmt/pull/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

PaulWessel commented 3 years ago

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

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

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

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

PaulWessel commented 3 years ago

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

PaulWessel commented 3 years ago

But makes no difference - same errors.

weiji14 commented 3 years ago

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

I probably need to install ipython then

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

PaulWessel commented 3 years ago

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

maxrjones commented 3 years ago

conda install ipython

PaulWessel commented 3 years ago

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

weiji14 commented 3 years ago

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

PaulWessel commented 3 years ago

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

maxrjones commented 3 years ago

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

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

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

Ok, PR https://github.com/GenericMappingTools/pygmt/pull/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 https://github.com/GenericMappingTools/pygmt/pull/750#discussion_r647908094).

maxrjones commented 2 years ago

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 https://github.com/GenericMappingTools/pygmt/pull/750#discussion_r647908094.

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

I will have a look at debugging this pygmt case.

maxrjones commented 2 years ago

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

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

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): ```python-traceback 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](https://user-images.githubusercontent.com/23487320/140242127-dbd728b1-97c6-4a21-9352-fa1fcf342176.png)

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

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: https://github.com/GenericMappingTools/gmt/blob/8c55a800d8103838fb94e70acf1b728cbc4740ff/src/gmt_api.c#L5242-L5327

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.

https://github.com/GenericMappingTools/gmt/blob/8c55a800d8103838fb94e70acf1b728cbc4740ff/src/gmt_grdio.c#L3235-L3264