Global earth relief data loaded as xarray doesn't work if projection center is not 0

GenericMappingTools / pygmt

A Python interface for the Generic Mapping Tools.

https://www.pygmt.org

BSD 3-Clause "New" or "Revised" License

747 stars 216 forks source link

Global earth relief data loaded as xarray doesn't work if projection center is not 0 #515

Closed seisman closed 4 years ago

seisman commented 4 years ago

Description of the problem

It may work in GMT 6.0.0 + PyGMT v0.1.2, or not. Currently, PyGMT master + GMT 6.1.0 gives me a wrong result.

Full code that generated the error

import pygmt
from pygmt.datasets import load_earth_relief

grid = load_earth_relief()
fig = pygmt.Figure()
fig.grdimage(grid, cmap="geo", region='d', projection="H10c")

fig.shift_origin(xshift="12c")
fig.grdimage(grid, cmap="geo", region='g', projection="H10c")

fig.shift_origin(xshift="-12c", yshift="6c")
fig.grdimage(grid, cmap="geo", region='d', projection="H120/10c")

# this one sometimes crashes
#fig.shift_origin(xshift="12c")
#fig.grdimage(grid, cmap="geo", region='g', projection="H120/10c")
fig.savefig("map.png")

Output

weiji14 commented 4 years ago

Not sure if #476 will fix this. I just refactored the code there but it's still giving the strange static-looking image in the gallery at https://pygmt-24935dva1.vercel.app/gallery/grid/track_sampling.html#sphx-glr-gallery-grid-track-sampling-py. This is with the pixel registered earth relief grid.

Sampling along tracks image distorted

Passing in the gridline registered makes it a bit better:

sphx_glr_track_sampling_001

seisman commented 4 years ago

@PaulWessel Please see if you have any idea why it fails. FYI, PyGMT is still using GMT_IN|GMT_IS_REFERENCE mode.

PaulWessel commented 4 years ago

Remind me: The xarray is being passed into GMT as a grid via a GMT_MATRIX structure?

seisman commented 4 years ago

Remind me: The xarray is being passed into GMT as a grid via a GMT_MATRIX structure?

Yes, GMT_IS_GRID|GMT_VIA_MATRIX.

PaulWessel commented 4 years ago

In gmtapi_import_grid:

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 */

Not completely reassuring... Could you try GMT_IN|GMT_IS_DUPLICATE and see if that works?

seisman commented 4 years ago

Using GMT_IN|GMT_IS_DUPLICATE, the GMT_Open_VirtualFile function returns a status code 32.

Using GMT_IN, it seems something wrong with the boundary condition:

PaulWessel commented 4 years ago

Sorry, what is different in settings between your TL panel here and the BL panel in the top post?

seisman commented 4 years ago

TL panel here: -Rd -JH120/10c and GMT_IN.

BL panel in the top post: -Rd -JH10c, and GMT_IN|GMT_IS_REFERENCE.

PaulWessel commented 4 years ago

OK, so -JH120/10.c does not to any rotation and messes up at -180 + 120 = 60W it seems.

seisman commented 4 years ago

The issue is caused by upstream GMT bugs, which was fixed in a series of upstream PRs (https://github.com/GenericMappingTools/gmt/pull/3813, https://github.com/GenericMappingTools/gmt/pull/3813, https://github.com/GenericMappingTools/gmt/pull/3829) and already merged into GMT's master and 6.1 branches.

GMT may release v6.1.1 this month. When GMT v6.1.1 is released, I think we can bump the minimum required GMT version to v6.1.1 and release PyGMT v0.2.0.

HELP NEEDED

PyGMT users who build GMT from source codes may try GMT's 6.1 branch (for the upcoming v6.1.1 release) and see if you can find any related bugs before GMT v6.1.1 is released.