GMTDataArrayAccessor: Support applying grid operations on the current xarray.DataArray object

[seisman]

This is a POC PR to show how GMTDataArrayAccessor can access GMT's grid-operation modules, i.e., grid.gmt.fill(**kwargs) instead of pygmt.grdfill(grid, **kwargs):

In [1]: import pygmt

In [2]: from pygmt.helpers.testing import load_static_earth_relief

In [3]: import numpy as np

In [4]: grid = load_static_earth_relief()

In [5]: grid[3:6, 3:5] = np.nan

In [6]: grid
Out[6]: 
<xarray.DataArray 'z' (lat: 14, lon: 8)> Size: 448B
array([[347.5, 344.5, 386. , 640.5, 617. , 579. , 646.5, 671. ],
       [383. , 284.5, 344.5, 394. , 491. , 556.5, 578.5, 618.5],
       [373. , 367.5, 349. , 352.5, 419.5, 428. , 570. , 667.5],
       [557. , 435. , 385.5,   nan,   nan, 496. , 519.5, 833.5],
       [561.5, 539. , 446.5,   nan,   nan, 553. , 726.5, 981. ],
       [310. , 521.5, 757. ,   nan,   nan, 524. , 686.5, 794. ],
       [521.5, 682.5, 796. , 886. , 571.5, 638.5, 739.5, 881.5],
       [308. , 595.5, 555.5, 556. , 580. , 770. , 927. , 920. ],
       [601. , 526.5, 535. , 299. , 398.5, 645. , 797.5, 964. ],
       [494.5, 488.5, 357. , 254.5, 286. , 484.5, 653.5, 930. ],
       [450.5, 395.5, 366. , 248. , 250. , 354.5, 550. , 797.5],
       [345.5, 320. , 335. , 292. , 207.5, 247. , 325. , 346.5],
       [349. , 313. , 325.5, 247. , 191. , 225. , 260. , 452.5],
       [347.5, 331.5, 309. , 282. , 190. , 208. , 299.5, 348. ]],
      dtype=float32)
Coordinates:
  * lon      (lon) float64 64B -54.5 -53.5 -52.5 -51.5 -50.5 -49.5 -48.5 -47.5
  * lat      (lat) float64 112B -23.5 -22.5 -21.5 -20.5 ... -12.5 -11.5 -10.5
Attributes:
    long_name:     elevation (m)
    actual_range:  [190. 981.]

In [7]: grid.gmt.fill(constantfill=20)
Out[7]: 
<xarray.DataArray 'z' (lat: 14, lon: 8)> Size: 448B
array([[347.5, 344.5, 386. , 640.5, 617. , 579. , 646.5, 671. ],
       [383. , 284.5, 344.5, 394. , 491. , 556.5, 578.5, 618.5],
       [373. , 367.5, 349. , 352.5, 419.5, 428. , 570. , 667.5],
       [557. , 435. , 385.5,  20. ,  20. , 496. , 519.5, 833.5],
       [561.5, 539. , 446.5,  20. ,  20. , 553. , 726.5, 981. ],
       [310. , 521.5, 757. ,  20. ,  20. , 524. , 686.5, 794. ],
       [521.5, 682.5, 796. , 886. , 571.5, 638.5, 739.5, 881.5],
       [308. , 595.5, 555.5, 556. , 580. , 770. , 927. , 920. ],
       [601. , 526.5, 535. , 299. , 398.5, 645. , 797.5, 964. ],
       [494.5, 488.5, 357. , 254.5, 286. , 484.5, 653.5, 930. ],
       [450.5, 395.5, 366. , 248. , 250. , 354.5, 550. , 797.5],
       [345.5, 320. , 335. , 292. , 207.5, 247. , 325. , 346.5],
       [349. , 313. , 325.5, 247. , 191. , 225. , 260. , 452.5],
       [347.5, 331.5, 309. , 282. , 190. , 208. , 299.5, 348. ]],
      dtype=float32)
Coordinates:
  * lat      (lat) float64 112B -23.5 -22.5 -21.5 -20.5 ... -12.5 -11.5 -10.5
  * lon      (lon) float64 64B -54.5 -53.5 -52.5 -51.5 -50.5 -49.5 -48.5 -47.5
Attributes:
    Conventions:   CF-1.7
    title:         
    history:       gmt grdfill @GMTAPI@-S-I-G-M-G-N-000000 -Ac20 -G@GMTAPI@-S...
    description:   
    actual_range:  [190. 981.]
    long_name:     z

Address #499 (comment).

Preview: https://pygmt-dev--3854.org.readthedocs.build/en/3854/api/generated/pygmt.GMTDataArrayAccessor.html

Grid-related modules

Processing (grid in, grid out)

• dimfilter
• grdclip
• grdcut
• grdfill
• grdfilter
• grdgradient
• grdhisteq: equalize_hist. The name comes from https://scikit-image.org/docs/stable/api/skimage.exposure.html#skimage.exposure.equalize_hist.
• grdproject
• grdsample
• grdtrack

Plotting

• grdcontour
• grdimage
• grdview

Misc

• grd2cpt
• grd2xyz
• grdinfo
• grdvolume

[weiji14]

Probably best to keep the method name as .gmt.grdclip instead of .gmt.clip, because there is also https://docs.generic-mapping-tools.org/6.5/clip.html (same with the other grd* methods below), and who knows if we might use pandas accessors at some point.

[seisman]

Even if we decide to implement the pandas accessors in the future (actually, I feel it's a good idea), there are still no conflicts since the two accessors are in different namespaces. E.g.,

• da.gmt.clip: Clip operation on a xr.DataArray grid/image
• df.gmt.clip: Clip operation on a pd.DataFrame object

[seisman]

I think this PR is ready for review. I guess I still need to add tests to increase code coverage?

[weiji14]

Shorter way to declare new methods is to put this in the __init__ method:

import functools

def __init__(self, xarray_obj: xr.DataArray):
    ...
    self.dimfilter = functools.partial(dimfilter, grid=self._obj)
    self.dimfilter.__doc__ = dimfilter.__doc__

This would preserve the full docs too, e.g. output from help(grid.gmt.dimfilter)

Help on partial in module functools:

functools.partial(<function dimfilter at 0x7fd07...e:  [190. 981.]
    long_name:     elevation (m))
    Directional filtering of grids in the space domain.
    
    Filter a grid in the space (or time) domain by
    dividing the given filter circle into the given number of sectors,
    applying one of the selected primary convolution or non-convolution
    filters to each sector, and choosing the final outcome according to the
    selected secondary filter. It computes distances using Cartesian or
    Spherical geometries. The output grid can optionally be generated as a
    subregion of the input and/or with a new increment using ``spacing``,
    which may add an "extra space" in the input data to prevent edge
    effects for the output grid. If the filter is low-pass, then the output
    may be less frequently sampled than the input. :func:`pygmt.dimfilter`
    will not produce a smooth output as other spatial filters
    do because it returns a minimum median out of *N* medians of *N*
    sectors. The output can be rough unless the input data are noise-free.
    Thus, an additional filtering (e.g., Gaussian via :func:`pygmt.grdfilter`)
    of the DiM-filtered data is generally recommended.
    
    Full option list at :gmt-docs:`dimfilter.html`
    
    **Aliases:**
    
    .. hlist::
       :columns: 3
    
       - D = distance
       - F = filter
       - I = spacing
       - N = sectors
       - R = region
       - V = verbose

There might be a nicer way to wrap things (maybe using https://docs.python.org/3/library/functools.html#functools.update_wrapper?), but haven't played around with it too much.

[seisman]

Tried your solution above and found two issues:

1. self.dimfilter is an attribute of the accessor, so it's not shown on the documentation (https://pygmt-dev--3854.org.readthedocs.build/en/3854/api/generated/pygmt.GMTDataArrayAccessor.html)
2. The help docs still show grid as its first parameter, which may be more confusing?

[seisman]

I've re-implemented this feature using functools.wraps (commit 397504b) and I think it works well.

• The methods are shown on the documentation: https://pygmt-dev--3854.org.readthedocs.build/en/3854/api/generated/pygmt.GMTDataArrayAccessor.html
• The help(grid.gmt.clip) docs is not perfect but is better:

Help on method grdfill in module pygmt.src.grdfill:

grdfill(outgrid: str | os.PathLike | None = None, constantfill: float | None = None, gridfill: str | os.PathLike | xarray.core.dataarray.DataArray | None = None, neighborfill: float | bool | None = None, splinefill: float | bool | None = None, inquire: bool = False, mode: str | None = None, *, hole=None, region=None, verbose=None, coltypes=None, **kwargs) -> xarray.core.dataarray.DataArray | numpy.ndarray | None method of pygmt.xarray.accessor.GMTDataArrayAccessor instance
    Interpolate across holes in a grid.

    Read a grid that presumably has unfilled holes that the user wants to fill in some
    fashion. Holes are identified by NaN values but this criteria can be changed via the
    ``hole`` parameter. There are several different algorithms that can be used to
    replace the hole values. If no holes are found the original unchanged grid is
    returned.

    Full GMT docs at :gmt-docs:`grdfill.html`.

    **Aliases:**

    .. hlist::
       :columns: 3

       - A = constantfill/gridfill/neighborfill/splinefill/mode
       - L = inquire
       - N = hole
       - R = region
       - V = verbose
       - f = coltypes

    Parameters
    ----------
    grid
        Name of the input grid file or the grid loaded as a
        :class:`xarray.DataArray` object.

        For reading a specific grid file format or applying basic data operations,
        see :gmt-docs:`gmt.html#grd-inout-full` for the available modifiers.

[Copilot]

Pull Request Overview

This PR introduces a new xarray accessor for grid operations on DataArray objects, simplifying the application of GMT’s grid-processing functionality directly on grids.

• Adds a set of accessor methods that wrap GMT grid-operation functions using a common _make_method utility.
• Adds test cases to validate the new grid operation methods (clip and equalize_hist) exposed via the accessor.

Reviewed Changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated no comments.

File	Description
pygmt/xarray/accessor.py	Introduces a helper method to wrap GMT grid functions and exposes grid operations as DataArray methods.
pygmt/tests/test_xarray_accessor.py	Adds tests to verify the new accessor methods work as expected.

Comments suppressed due to low confidence (2)

pygmt/xarray/accessor.py:259

• [nitpick] Although the 'filter' method is appropriately scoped to the accessor, its name coincides with the built-in Python function 'filter'. Consider renaming it (e.g., 'apply_filter') to avoid any potential confusion.

    filter = _make_method(grdfilter)

pygmt/tests/test_xarray_accessor.py:218

• The comment references 'test_grdtest_grdclip_no_outgrid', which appears to be a typo. Please update the reference to clearly indicate the correct test name or remove the ambiguous reference for clarity.

    This test is adapted from the `test_grdtest_grdclip_no_outgrid` test.