Cache gmt accessor info for non-virtualfile grid/image inputs

[weiji14]

Description of proposed changes

Ensure that using load_dataarray/open_dataarray(..., engine="gmt", ...) caches the GMTDataArrayAccessor info, so that the correct registration/gtype is returned even if the original source NetCDF file is deleted.

TODO:

• Add regression test
• Fix bug by caching the accessor values properly

Fixes #4005

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If wrapping a new module, open a 'Wrap new GMT module' issue and submit reasonably-sized PRs.
• If adding new functionality, add an example to docstrings or tutorials.

Slash Commands

You can write slash commands (/command) in the first line of a comment to perform
specific operations. Supported slash command is:

• /format: automatically format and lint the code

[weiji14]

Originally posted by @seisman in #3932 (comment)

Actually, it's likely that the accessor information will be lost when converting via to_dataset.

So we can set the accessor info at the xr.Dataset level, but not sure if it's a good idea because then xr.Dataset would have GMT-functions enabled...

[seisman]

The test is still failing. I guess that the accessor information is lost in both dataarray->dataset and dataset->dataarray conversions.

When calling xr.open_dataarray with the "gmt" engine, it actually calls the GMTBackendEntry's open_dataset method first, so the conversion from dataset to dataarray cannot be avoided.

[seisman]

>>> import pygmt
>>> import xarray as xr

>>> grid = xr.open_dataarray("@static_earth_relief.nc", engine="gmt", raster_kind="grid")

>>> grid.encoding["source"]
np.str_('/home/seisman/.gmt/cache/static_earth_relief.nc')

>>> grid.gmt.registration, grid.gmt.gtype
(<GridRegistration.PIXEL: 1>, <GridType.GEOGRAPHIC: 1>)

>>> # Set source encoding to a fake path, similar to a temporary file being deleted
>>> grid.encoding["source"] = "/path/to/nothing"  
>>> grid.gmt.registration, grid.gmt.gtype
(<GridRegistration.PIXEL: 1>, <GridType.GEOGRAPHIC: 1>)

>>> # Convert a dataarray to a dataset. The source encoding is not kept.
>>> da = grid.to_dataset()
>>> da.encoding
{}

>>> # Convert a dataset to a dataarray.
>>> # Xref: https://github.com/pydata/xarray/blob/fb49a3b0a83272467aa01ac189b95e9e7ebe01aa/xarray/backends/api.py#L950
>>> (ds, ) = da.data_vars.values()
>>> ds.encoding
{'source': '/path/to/nothing'}

>>> ds.gmt.registration, ds.gmt.gtype
(<GridRegistration.GRIDLINE: 0>, <GridType.CARTESIAN: 0>)

[seisman]

As shown above, the accessor information is lost during the dataset->dataarray conversion, and we rely on the source encoding to decide the registration and gtype, which default to 0 for temporary files being deleted.

I don't think we have any solutions here, and need to document the limitation and recommend:

1. Avoiding using temporary files if possible
2. If temporary files cannot be avoided, use _ = grid.gmt before deleting the temporary file

[weiji14]

Yeah, I found this upstream issue - pydata/xarray#3268 that accessors are not meant to cache a state. I do wish though that the @xr.register_dataarray_accessor decorator could work by us returning an xr.DataArray directly (instead of having to return an xr.Dataset that then gets converted to an xr.DataArray).

Agree that we should document that temporary files should be avoided and/or trigger caching using _ = grid.gmt, then close #4005 as wontfix (with a workaround).

[weiji14]

I do wish though that the @xr.register_dataarray_accessor decorator could work by us returning an xr.DataArray directly (instead of having to return an xr.Dataset that then gets converted to an xr.DataArray).

Opened an upstream issue on xarray for this -> pydata/xarray#10562