Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
eelcodoornbos opened this issue Nov 19, 2022 · 6 comments
Labels
feature request New feature wanted

Comments

@eelcodoornbos
Copy link

eelcodoornbos commented Nov 19, 2022

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: #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.

@eelcodoornbos eelcodoornbos added the feature request New feature wanted label Nov 19, 2022
@welcome
Copy link

welcome bot commented Nov 19, 2022

👋 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
Copy link
Member

seisman commented Nov 20, 2022

@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
Copy link
Member

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

@weiji14
Copy link
Member

weiji14 commented Nov 21, 2022

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

pygmt/pygmt/accessors.py

Lines 28 to 39 in d8c816a

def __init__(self, xarray_obj):
self._obj = xarray_obj
try:
self._source = self._obj.encoding["source"] # filepath to NetCDF source
# Get grid registration and grid type from the last two columns of
# the shortened summary information of `grdinfo`.
self._registration, self._gtype = map(
int, grdinfo(self._source, per_column="n").split()[-2:]
)
except (KeyError, ValueError):
self._registration = 0 # Default to Gridline registration
self._gtype = 0 # Default to Cartesian grid type

Currently we do try to use grdinfo to determine the grid registration and type if the xarray.DataArray was loaded from a NetCDF (see #500 (comment)). 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
Copy link
Member

seisman commented Dec 31, 2022

@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
Copy link
Member

Yes, I would think that is reasonable

@seisman seisman changed the title grdimage to infer geographic grid when provided with xarray input Infer grid registration (gridline vs pixel) and grid type (geographic or Cartesian) for xarray.DataArray input Feb 18, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
feature request New feature wanted
Projects
None yet
Development

No branches or pull requests

4 participants