PyGMT: Gap in grdimage of xarray data at zero longitude

eelcodoornbos · November 11, 2022, 6:29pm

I’m trying to plot an xarray grid on a map with pygmt. Here’s a constructed example using the same array data structure plotted with pygmt as well as saved as netcdf and plotted using GMT in the OS shell.

import os
import numpy as np
import xarray as xr
import pygmt
from IPython.display import Image

latitude = np.arange(-90,90,2, dtype='float')
longitude = np.arange(0,360,4, dtype='float')
xx, yy = np.meshgrid(longitude, latitude)
fakedata = np.cos(2*np.deg2rad(xx)) * np.cos(4*np.deg2rad(yy))

xrdata = xr.DataArray(data=fakedata,
                      coords={'longitude': longitude,
                              'latitude': latitude},
                      dims=('latitude', 'longitude'))

fig = pygmt.Figure()
fig.grdimage(xrdata, projection='G0/30/15c', region='d', dpi=300, interpolation='n')
fig.savefig('fakedata_pygmt.png')

xrdata.to_netcdf('fakedata.nc')
os.system("gmt grdimage fakedata.nc -JG0/30/15c -Rg -E300 -nn > fakedata_shell.ps;"+
          "gmt psconvert -Tg -A fakedata_shell.ps")

The output is available here:

The pygmt version contains a gap at 0 zero longitude. I would have liked to use pygmt to create a plot that wraps in longitude like in the OS shell version. Is this possible? Is this a bug in pygmt, or do I need to do something to the input in pygmt to achieve this?

yvonnefroehlich · November 11, 2022, 7:17pm

Hello @eelcodoornbos,

please note, that numpy.arange excludes the upper boundary, i.e.

longitude = np.arange(0,360,4, dtype="float")

stops at 356, not at 360. This is probably the reason for your gap at 360° or 0° longitude.

And similar for the gap at latitude 90° North, as

latitude = np.arange(-90,90,2, dtype="float")

stops at 88, not at 90.

P. S.
You can format your script as code by placing three backticks (```) in the line before and after the block with the code:

```
# your script formated as code
```

eelcodoornbos · November 11, 2022, 7:31pm

Thank you Yvonne, yes I knew of the np.arange behaviour, but in my example I’m recreating the grid of the output from a global atmospheric/ionospheric model, so if that is the reason, is the model in error concerning the longitude output? Does this mean I have to repeat the data at 0 and 360 longitude in input to PyGMT? Is this a documented convention? The incomplete latitude grid is my mistake in this example, 90N is in the model output.

And why would this behaviour be different between command-line GMT and pyGMT? I think the command-line GMT behaviour makes more sense. Grdimage uses interpolation (nearest-neighbour in this case), and I expected GMT to know that for a pixel at ~359 deg longitude, the nearest neighbour is at 0 deg longitude in this grid. It seems to do so in the command-line version.

eelcodoornbos · November 11, 2022, 8:19pm

When I change the code like this:

latitude = np.arange(-90,91,2, dtype='float')
longitude = np.arange(0,361,4, dtype='float')

The result indeed looks better in pyGMT, but it gives this warning, which I don’t understand:

grdimage [WARNING]: Longitude range too small; geographic boundary condition changed to natural.

yvonnefroehlich · November 12, 2022, 11:49am

There seems to be a difference in PyGMT between passing a xr.DataArray or a netCDF file through the grid parameter of pygmt.grdimage. When using the netCDF file the gap at 360° or 0° longitude does not occur. Currently, I do not know why this is.

import os
import numpy as np
import xarray as xr
import pygmt


latitude = np.arange(-90,91,2, dtype="float")
longitude = np.arange(0,360,4, dtype="float")

xx, yy = np.meshgrid(longitude, latitude)
fakedata = np.cos(2*np.deg2rad(xx)) * np.cos(4*np.deg2rad(yy))

xrdata = xr.DataArray(
    data=fakedata,
    coords={"longitude":longitude, "latitude":latitude},
    dims=("latitude", "longitude")
)
xrdata.to_netcdf("fakedata.nc")

# -----------------------------------------------------------------------------
# use xarray DataArray as input
fig = pygmt.Figure()

fig.grdimage(
    grid=xrdata,
    projection="G0/30/15c",
    region="d",
    dpi=300,
    interpolation="n",
)

fig.show()
# fig.savefig(fname="fakedata_pygmt.png")

# -----------------------------------------------------------------------------
# use netcdf file as input
fig = pygmt.Figure()

fig.grdimage(
    grid="fakedata.nc",
    projection="G0/30/15c",
    region="d",
    dpi=300,
    interpolation="n",
)

fig.show()
# fig.savefig(fname="fakedata_pygmt_nc.png")

output figure when using the xr.DataArray:

output figure when using the netcdf file:

P. S.: Thanks for formating your code @eelcodoornbos

eelcodoornbos · November 12, 2022, 1:26pm

That’s a good workaround using NetCDF, but of course it should work with XArray input as well. There seems to be a more general problem with grdimage and the perspective projection in PyGMT using XArray input instead of a NetCDF file. I’ve now tried a geostationary satellite perspective:
projection="G-47/0/35792.29/0/0/0/50/50/10c"
This results in the entire western hemisphere becoming empty when using XArray, not when using the NetCDF file saved via XArray. When changing the central latitude of the projection (the first zero above) to about 9 or higher, the issue goes away, perhaps because the pole comes in view.

I wonder if internally pygmt does some sort of preprocessing of the longitude coordinates in XArray grid data that it does not apply to data read from a NetCDF file (or vice versa). If so, I expect that there’s a bug in that part of the code. Perhaps that code is also where the warning comes from. I have not been able to find it yet.

yvonnefroehlich · November 12, 2022, 6:53pm

Ping @seisman, @weiji14, and @maxrjones for thoughts and help on this, please.
Maybe it makes sees to open a related issue in the GitHub repository at GenericMappingTools/pygmt: A Python interface for the Generic Mapping Tools. (github.com)?

seisman · November 13, 2022, 7:29am

xrdata = xr.DataArray(data=fakedata,
                      coords={'longitude': longitude,
                              'latitude': latitude},
                      dims=('latitude', 'longitude'))

The main reason is that, PyGMT has no idea if xrdata is a Cartesian grid or a geographic grid. So the solution is telling PyGMT that this grid is a geographic grid by setting the GMT’s xarray accessor (https://www.pygmt.org/dev/api/generated/pygmt.GMTDataArrayAccessor.html#pygmt.GMTDataArrayAccessor):

xrdata = xr.DataArray(data=fakedata,
                      coords={'longitude': longitude,
                              'latitude': latitude},
                      dims=('latitude', 'longitude'))
xrdata.gmt.gtype = 1

Here gtype = 0 means Cartesian grid and gtype = 1 means geographic grid.

eelcodoornbos · November 14, 2022, 8:26am

Thanks @seisman, I had been looking for something like that. I did some more testing and the gap in the original problem only appears when both registration and gtype are set to 0. In all other cases, there is no gap. This makes sense.

When plotting the NetCDF file, there is never a gap, however. When using grdinfo and playing around with variants of xarray data saved to NetCDF files, I can also see that GMT can infer whether a geographic or Cartesian grid type is used, by whether or not there are coordinates named “longitude” and “latitude” (COARDS conventions I guess).
I think it would be helpful if grdimage would infer this as well, when fed with a fresh xarray structure, because the current solution of manually setting the gtype is not very easy to discover. I’m a 20+ year GMT user trying to convince some of my colleagues to try out pygmt, but these are the type of things that make them stop using it.

At the same time, the problem with grdimage using an xarray in the low-latitude perspective (geostationary satellite) projection seems to be a separate problem, which only gets worse when setting the gtype. When gtype=0 (Cartesian) only the Eastern hemisphere is plotted (as was shown above). When I now set gtype=1 (geographic), it does not even seem to read the z-values correctly, as it triggers a GMTCLibError with the following message:

grdimage [ERROR]: Failed to read CPT (null).
psconvert [ERROR]: System call [gs -q -dNOSAFER -dNOPAUSE -dBATCH -sDEVICE=bbox -DPSL_no_pagefill -dMaxBitmap=2147483647 -dUseFastColor=true '/Users/doornbos/.gmt/sessions/gmt_session.43073/gmt_112.ps-' 2> '/Users/doornbos/.gmt/sessions/gmt_session.43073/psconvert_43073c.bb'] returned error 256.
Error: /undefined in PSL_movie_label_completion
Operand stack:
   PSL_eoc
Execution stack:
   %interp_exit   .runexec2   --nostringval--   --nostringval--   --nostringval--   2   %stopped_push   --nostringval--   --nostringval--   --nostringval--   false   1   %stopped_push   1990   1   3   %oparray_pop   1989   1   3   %oparray_pop   1977   1   3   %oparray_pop   1833   1   3   %oparray_pop   --nostringval--   %errorexec_pop   .runexec2   --nostringval--   --nostringval--   --nostringval--   2   %stopped_push   --nostringval--
Dictionary stack:
   --dict:770/1123(ro)(G)--   --dict:0/20(G)--   --dict:75/200(L)--   --dict:147/250(L)--
Current allocation mode is local
Last OS error: Invalid argument
Current file position is 20517
Error plotting reg=0,gtype=1

I find it weird that this error depends on the projection central latitude. When explicitly passing a CPT with the makecpt command, this error message goes away, but it results in an empty plot. Only the grdimage incantation which provides a NetCDF-filename as input seems to work with these projection parameters.

seisman · November 14, 2022, 11:40pm

When plotting the NetCDF file, there is never a gap, however. When using grdinfo and playing around with variants of xarray data saved to NetCDF files, I can also see that GMT can infer whether a geographic or Cartesian grid type is used, by whether or not there are coordinates named “longitude” and “latitude” (COARDS conventions I guess).

Yes, GMT is more clever about guessing the grid type and registration, while PyGMT doesn’t have the ability now. c

I think it would be helpful if grdimage would infer this as well, when fed with a fresh xarray structure, because the current solution of manually setting the gtype is not very easy to discover. I’m a 20+ year GMT user trying to convince some of my colleagues to try out pygmt, but these are the type of things that make them stop using it.

Could you please open a feature request in the PyGMT repository?

At the same time, the problem with grdimage using an xarray in the low-latitude perspective (geostationary satellite) projection seems to be a separate problem, which only gets worse when setting the gtype. When gtype=0 (Cartesian) only the Eastern hemisphere is plotted (as was shown above). When I now set gtype=1 (geographic), it does not even seem to read the z-values correctly, as it triggers a GMTCLibError with the following message:

Yes, there are a series of related issues (https://github.com/GenericMappingTools/pygmt/issues/732, Global earth relief data loaded as xarray doesn't work if projection center is not 0 · Issue #515 · GenericMappingTools/pygmt · GitHub, Cylindrical projections 'Q' plot incorrectly in pygmt when using xarrays · Issue #390 · GenericMappingTools/pygmt · GitHub, Geographic xarray.DataArray grid is plotted as Cartesian grid · Issue #1172 · GenericMappingTools/pygmt · GitHub, grdimage does not correctly plot xarray DataArrays when they don't contain the redundant longitude at 360 E · Issue #375 · GenericMappingTools/pygmt · GitHub). Some are also fixed while some are not. Could you please open a bug report if your issue is not covered by the previously reported issues?