grdimage does not plot subsets of data correctly on a global map

[MarkWieczorek]

Problem

I have a subset of data that I would like to plot on a global projection (such as mollweide). This works with GMT, but pygmt provides incorrect output.

Example

Let's create a global dataset of random numbers on a grid with 1 degree spacing in latitude and longitude:

lats = np.arange(90., -91., -1.)
lons = np.arange(0., 361., 1.)
data = np.random.rand(181, 361)
da = xr.DataArray(data, coords=[('latitude', lats), ('longitude', lons)])

Now let's plot only the region from 90N-30S and 0E-220E with a Mollweide projection centered on 180 E

figure = pygmt.Figure()
pygmt.makecpt(series=[0., 1., 0.1], cmap='vik', continuous=True)
figure.grdimage(da[0:121, 0:221], region='g', projection='W180/8i', frame=True)
figure.show(method='external')

This is as expected, but now do the same, but centered on 0 E

figure = pygmt.Figure()
pygmt.makecpt(series=[0., 1., 0.1], cmap='vik', continuous=True)
figure.grdimage(da[0:121, 0:221], region='g', projection='W0/8i', frame=True)
figure.show(method='external')

And as a second example, here I want to plot the global data, but using a dataset the doesn't contain the redundant data at 360 E

figure = pygmt.Figure()
pygmt.makecpt(series=[0., 1., 0.1], cmap='vik', continuous=True)
figure.grdimage(da[0:181, 0:360], region='g', projection='W0/8i', frame=True)
figure.show(method='external')

This is a complete failure, but changing the central meridian to 180 degrees works!

figure = pygmt.Figure()
pygmt.makecpt(series=[0., 1., 0.1], cmap='vik', continuous=True)
figure.grdimage(da[0:181, 0:360], region='g', projection='W180/8i', frame=True)
figure.show(method='external')

[seisman]

@PaulWessel I think this is also an upstream GMT bug, possibly related to GenericMappingTools/gmt#4335. Please see if you can help.

To reproduce the issue using the above codes, you need to import the following packages first:

import pygmt
import numpy as np
import xarray as xr

[seisman]

Ping @PaulWessel

[weiji14]

Still an issue on GMT 6.2.0rc1

import numpy as np
import pygmt
import xarray as xr

lats = np.arange(90.0, -91.0, -1.0)
lons = np.arange(0.0, 361.0, 1.0)
data = np.random.rand(181, 361)
da = xr.DataArray(data, coords=[("latitude", lats), ("longitude", lons)])

fig = pygmt.Figure()
pygmt.makecpt(series=[0.0, 1.0, 0.1], cmap="vik", continuous=True)
fig.grdimage(da[0:121, 0:221], region="g", projection="W0/8i", frame=True)
fig.savefig("temp.png")
fig.show()

[PaulWessel]

Busy with taxes this weekend but will get on this - sorry it slipped off the radar. Seems like I fixed a very related problem not that long ago, no @seisman ?

[PaulWessel]

Not sure if https://docs.generic-mapping-tools.org/dev/debug.html#debug-pygmt-in-xcode-on-macos has all required steps. After git pull, make install, conda activate pygmt I try (from pygmt directory)

python ~/bad.py (this is the script posted above)
Traceback (most recent call last):
  File "/Users/pwessel/bad.py", line 2, in <module>
    import pygmt
  File "/Users/pwessel/GMTdev/pygmt/pygmt/__init__.py", line 28, in <module>
    from pygmt.figure import Figure, set_display
  File "/Users/pwessel/GMTdev/pygmt/pygmt/figure.py", line 9, in <module>
    import IPython
ModuleNotFoundError: No module named 'IPython'

[seisman]

Not sure if docs.generic-mapping-tools.org/dev/debug.html#debug-pygmt-in-xcode-on-macos has all required steps. After git pull, make install, conda activate pygmt I try (from pygmt directory)

python ~/bad.py (this is the script posted above)
Traceback (most recent call last):
  File "/Users/pwessel/bad.py", line 2, in <module>
    import pygmt
  File "/Users/pwessel/GMTdev/pygmt/pygmt/__init__.py", line 28, in <module>
    from pygmt.figure import Figure, set_display
  File "/Users/pwessel/GMTdev/pygmt/pygmt/figure.py", line 9, in <module>
    import IPython
ModuleNotFoundError: No module named 'IPython'

I think you found a PyGMT bug which was introduced recently.

Your temporary fix is:

conda activate python
conda install ipython

[PaulWessel]

So I am on miniconda:

(pygmt) pwessel@macnut:~/GMTdev/pygmt-> conda activate python
Could not find conda environment: python
You can list all discoverable environments with conda info --envs.

[seisman]

Sorry, typo. I meant:

conda activate pygmt
conda install ipython

[PaulWessel]

OK, that installed a bunch of stuff including ipython. But still some PATH issues perhaps (remember, I do this so infrequently I may not have things set up right):

Traceback (most recent call last):
File "/Users/pwessel/bad.py", line 2, in
import pygmt
File "/Users/pwessel/GMTdev/pygmt/pygmt/init.py", line 28, in
from pygmt.figure import Figure, set_display
File "/Users/pwessel/GMTdev/pygmt/pygmt/figure.py", line 9, in
import IPython
ModuleNotFoundError: No module named 'IPython'

[seisman]

What's the output of which python? Just want to make sure that you're running the one in the pygmt environment.

[PaulWessel]

/Users/pwessel/opt/miniconda3/bin/python

[seisman]

The path is incorrect. After activating the pygmt environment, you should be using /Users/pwessel/opt/miniconda3/envs/pygmt/bin/python. Did you forget to run conda activate pygmt?

[PaulWessel]

No, and just run it again. Same path coming back.

[seisman]

It doesn't make sense to me.

At line 10 of file /Users/pwessel/GMTdev/pygmt/pygmt/figure.py, changing KeyError to ModuleNotFoundError, then try to run the bad.py script, does it work for you?

[PaulWessel]

That runs fine. I get a ton of debug memory leak messages that should be fixed in general but the script works and produces the failed image.

[weiji14]

Ok, some testing on GMT 6.2.0rc2 (install using conda install -c conda-forge/label/dev gmt=6.2.0rc2):

# Random data setup
import numpy as np
import pygmt
import xarray as xr

lats = np.arange(90.0, -91.0, -1.0)
lons = np.arange(0.0, 361.0, 1.0)
data = np.random.rand(181, 361)
da = xr.DataArray(data, coords=[("latitude", lats), ("longitude", lons)])

Case 1

Previous whitespaces on the left and right are now coloured in correctly.

fig = pygmt.Figure()
pygmt.makecpt(series=[0.0, 1.0, 0.1], cmap="vik", continuous=True)
fig.grdimage(da[0:121, 0:221], region="g", projection="W0/8i", frame=True)
fig.savefig("temp1.png")
fig.show()

Case 2

This previously failed with a complete blue image. Now the random data shows up, but there is still a white line running through the middle

figure = pygmt.Figure()
pygmt.makecpt(series=[0., 1., 0.1], cmap='vik', continuous=True)
figure.grdimage(da[0:181, 0:360], region='g', projection='W0/8i', frame=True)
figure.savefig("temp2.png")
figure.show()

[PaulWessel]

You go 0 to 181 in latitude but not 0 to 361 in longitude? Aren't you missing a column?

[weiji14]

You go 0 to 181 in latitude but not 0 to 361 in longitude? Aren't you missing a column?

Ah right, thanks Paul. So we do need a redundant column to wrap around the dateline. Using figure.grdimage(da[0:181, 0:361], region="g", projection="W0/8i", frame=True) removes the middle white line (though it complains with grdimage [WARNING]: Longitude range too small; geographic boundary condition changed to natural.):

Ok to close this issue @seisman since it's resolved with GMT 6.2.0rc2? Or do we want to do some more work/testing on this?

[seisman]

Ok to close this issue @seisman since it's resolved with GMT 6.2.0rc2? Or do we want to do some more work/testing on this?

I think we need to add a test for it although it's an upstream bug.

[weiji14]

Just a note that doing a slice of an xarray.DataArray grid resets the coordinate system from Geographic type to Cartesian type as mentioned in #524 (comment) and this may be part of the issue (specifically for the 'white line in the middle' case 2 in #732 (comment)).