GMTDataArrayAccessor: Enable grid operations on the current xarray.DataArray object directly

[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(mode="c20")
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?