Infer grid registration (gridline vs pixel) and grid type (geographic or Cartesian) for xarray.DataArray input

[eelcodoornbos]

Currently, grdimage in pygmt can accept a NetCDF file name or an xarray structure as input for the grid to be plotted. When providing a NetCDF file, the code apparently checks whether the grid is geographic (has latitude/longitude dimensions), and behaves accordingly, by closing the grid at 0/360 deg longitude. When providing an xarray structure, however, such a check is not performed, and the user has to set the .gmt.gtype property on the xarray data to 1 to make pygmt plot the data as a geographic grid.

The topic was discussed on the forum, where a code example is provided: https://forum.generic-mapping-tools.org/t/pygmt-gap-in-grdimage-of-xarray-data-at-zero-longitude/3431/9

I've also just found this closely related (closed) issue: https://github.com/GenericMappingTools/pygmt/issues/375

A temporary solution might be to improve the documentation on how the user should set the xarray .gmt.gtype property in the documentation, perhaps by including an example in the Tuturial related to grdimage (https://www.pygmt.org/latest/tutorials/advanced/earth_relief.html#sphx-glr-tutorials-advanced-earth-relief-py), or by adding a separate tutorial on how NetCDF files and xarray structures can be created from Python and are handled together with PyGMT.

Are you willing to help implement and maintain this feature? No, I'm not at the moment able to invest the necessary time in working out the inner workings of the code, especially since this seems related to a complex long-running series of issue. I could help with the documentation, though.

[welcome[bot]]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. You might also want to take a look at our contributing guidelines and code of conduct.

[seisman]

@PaulWessel I remember that GMT has a function or lines of code that checks if a grid is geographic or Cartesian, but I can't find where it is. Could you please point me to the function if you remember?

[PaulWessel]

In usual situations (a grid has ben read or -R -J processed) then we use the macro _gmt_M_isgeographic (GMT, GMT_IN) which returns true of we have a geographic situation. There is also _gmt_M_iscartesian (GMT, GMT_IN) as well.

[weiji14]

Hmm, depending on how that geographic/cartesian logic works, we could probably add it to the logic around here:

https://github.com/GenericMappingTools/pygmt/blob/d8c816a6065f151003eaf021d83b0359c9995108/pygmt/accessors.py#L28-L39

Currently we do try to use grdinfo to determine the grid registration and type if the xarray.DataArray was loaded from a NetCDF (see https://github.com/GenericMappingTools/pygmt/pull/500#discussion_r453117559). But for a newly constructed xarray.DataArray, the default is gridline registration and Cartesian. We discussed in #487 on whether the default should change from gridline to pixel registration, and from Cartesian to Geographic.

If I'm hearing correctly, resolving this issue would require

1. adding some smarter logic around determining whether the xarray.DataArray is in fact Geographic or Cartesian, rather than just falling back to the default
2. Documenting the .gtype/.registration flag in a tutorial to be clear for xarray users that might not realize what the GMT default is

[seisman]

@PaulWessel I remember that GMT has a function or lines of code that checks if a grid is geographic or Cartesian, but I can't find where it is. Could you please point me to the function if you remember?

I just found this function gmtnc_grd_info, and found a long explanation about how GMT determines the grid registration for reading a netCDF file that doesn't have the node_offset attributes. I feel the same logic can also be applied to a user-created xarray.DataArray which usuallay only contains x/y/z arrays.

  /* Explanation for the logic below: Not all netCDF grids are proper COARDS grids and hence we sometime must guess
   * regarding the settings.  The x and y coordinates may be written as arrays, which reflect the positions of the
   * nodes.  There may also be the attributes actual_range which specifies the x range.  These will differ if the
   * grid is pixel registered, hence when both are present we use this difference to detect a pixel grid. However,
   * some external tools such as xarray may slice a grid but not update the attributes.  In this case the actual_range
   * may have an initial range that is no longer the case.  We have added a check if these differ by more than a
   * half grid increment.  If not then we can trust it.  If actual_range is missing then we have to guess the registration
   * which we do by checking if the start coordinate is an integer multiple of the increment.  If not, we guess pixel registration
   * but cannot know if this is the case unless the adjusted coordinates in x has a range of 360 and in y a range of 180.
   * Finally, if there is no array just the actual_range, then we cannot tell the registration from the range but try
   * and leave it as gridline registration. */

I also found the gmtlib_get_grdtype function, which determines if a grid is Cartesian or Geographic based on x and y ranges.

It seems that we can apply the same logic to xarray.DataArray data.

[PaulWessel]

Yes, I would think that is reasonable