Cylindrical projections 'Q' plot incorrectly in pygmt when using xarrays

[MarkWieczorek]

Pygmt incorrectly plots xarrays when using cylindrical projections ('Q').

This is not an issue with gmt: if the dataarray is first exported to a netcdf file, pygmt works as intended when using the file. Thus, this appears to be a problem with how pygmt treats internal xarrays, and only for the case of cylindrical projections (the most basic of all the projections!). This is demonstrated in the following script.

The script generates a simple dataset, and the plots the following:

• upper row: 'Q0/0/6i' and 'Q180/0/6i' reading data from a file.
• middle row: 'W0/6i' and 'W180/6i' using internal xarray.
• bottom row: 'Q0/0/6i' and 'Q180/0/6i' using internal xarray.

As you see, the upper two rows plot the data correctly. However, for the bottom left image, the data are not plotted with the correct central_longitude. Furthermore, there is a problem with the first longitude band not being plotted.

It is possible that this is related to this issue: https://github.com/GenericMappingTools/pygmt/issues/375

# Creata a data array in gridline coordinates of sin(lon) * cos(lat)
interval = 10
lat = np.arange(90, -90-interval, -interval)
lon = np.arange(0, 360+interval, interval)
longrid, latgrid = np.meshgrid(lon, lat)
data = np.sin(np.deg2rad(longrid)) * np.cos(np.deg2rad(latgrid))

# create a DataArray
dataarray = xr.DataArray(data, coords=[('latitude', lat,
                                       {'units': 'degrees_north'}),
                                       ('longitude', lon,
                                       {'units': 'degrees_east'})], 
                         attrs = {'actual_range': [-1, 1]})
dataset = dataarray.to_dataset(name='dataarray')
dataset.to_netcdf('test.grd')

# create projected images using a cylindrical and mollweide projection.
fig = pygmt.Figure()
fig.grdimage(dataset.dataarray, region='g', projection='Q0/0/6i', frame='a90f30g30')
fig.shift_origin(xshift='7i')
fig.grdimage(dataset.dataarray, region='g', projection='Q180/0/6i', frame='a90f30g30')
fig.shift_origin(xshift='-7i', yshift='4i')
fig.grdimage(dataset.dataarray, region='g', projection='W0/6i', frame='a90f30g30')
fig.shift_origin(xshift='7i')
fig.grdimage(dataset.dataarray, region='g', projection='W180/6i', frame='a90f30g30')
fig.shift_origin(xshift='-7i', yshift='4i')
fig.grdimage('test.grd', region='g', projection='Q0/0/6i', frame='a90f30g30')
fig.shift_origin(xshift='7i')
fig.grdimage('test.grd', region='g', projection='Q180/0/6i', frame='a90f30g30')

fig.savefig('test.png')
fig.show()

[MarkWieczorek]

This is just an update about the problem pygmt is having with xarray grids. I just tried the above example with the latest commit on master, and am now getting the following error:

In [23]: fig = pygmt.Figure()
    ...: fig.grdimage(dataset.dataarray, region='g', projection='Q0/0/6i', frame='a90f3
    ...: 0g30')
grdimage [ERROR]: Passing zmax <= zmin prevents automatic CPT generation!
grdimage [ERROR]: Failed to read CPT (null).

[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.

@MarkWieczorek If you can build GMT from source codes, please try the GMT's 6.1 branch and see if you still have any grid related issues.

[MarkWieczorek]

After reinstalling gmt via brew from HEAD (brew install gmt --HEAD), this now appears to work for me too. Thanks!

[seisman]

I'm reopening the issue so that we won't forget to add some tests for this issue.

[weiji14]

I'm reopening the issue so that we won't forget to add some tests for this issue.

Just for reference (since it took me a while to remember), we did add some tests for Q-type Cylindrical Equidistant projections already at https://github.com/GenericMappingTools/pygmt/pull/560#discussion_r484011483, but certain lon0/lat0 combinations (e.g. Q123/0, Q123/30, Q180/0, Q180/30) still don't work (as of GMT 6.2.0rc2). We'll close this issue down once those Q combinations work in PyGMT.

[weiji14]

Yet another half-yearly update. Upstream GMT issues to track this are at https://github.com/GenericMappingTools/gmt/issues/4335#issuecomment-709729148 and https://github.com/GenericMappingTools/gmt/issues/4440#issuecomment-724961749. Might be a longstanding issue with pixel registered grids?

[MarkWieczorek]

Just checked in after a long hiatus from coding, and it seems that this still isn't working (unfortunately). I suspect that this is unrelated to the GMT thread that you noted above, as the problem is simply related to setting the prime meridian.

[MarkWieczorek]

For info, I just checked that this problem still persists with pygmt 0.12. I was hoping that the virtual file changes would have fixed this :(

[weiji14]

Yes, the virtualfile changes in v0.12.0 were only for the output (i.e. going from GMT virtualfile -> xarray). To fix this, we'll need to handle the reverse direction (i.e. going from xarray -> GMT virtualfiles). Keep an eye on #3099 if you're interested :wink: