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

[seisman]

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]

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.

Passing in the gridline registered makes it a bit better:

[seisman]

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

[PaulWessel]

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

[seisman]

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]

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]

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]

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

[seisman]

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

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

[PaulWessel]

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

[seisman]

The issue is caused by upstream GMT bugs, which was fixed in a series of upstream PRs (GenericMappingTools/gmt#3813, GenericMappingTools/gmt#3813, GenericMappingTools/gmt#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.