Passing xarray.DataArray grids through data processing functions should preserve dimension names

[weiji14]

Description of the desired feature

Currently, the dimension names of the xarray.DataArray grid outputs appear to be hardcoded to 'x', 'y' and 'z', instead of using the original dimension names of the input xarray.DataArray. To illustrate

grid = xr.DataArray(
    data=[
        [0, 1, 2],
        [3, 4, 5],
        [6, 7, 8],
    ],
    coords=dict(height=[1.5, 2.5, 3.5], width=[4.5, 5.5, 6.5]),
    dims=["height", "width"],
)
print(grid)
# <xarray.DataArray (height: 3, width: 3)>
# array([[0, 1, 2],
#        [3, 4, 5],
#        [6, 7, 8]])
# Coordinates:
#   * height   (height) float64 1.5 2.5 3.5
#   * width    (width) float64 4.5 5.5 6.5

Run it through a data processing function like grdhisteq.equalize_grid

newgrid = pygmt.grdhisteq.equalize_grid(grid=grid)
print(newgrid)
# <xarray.DataArray 'z' (y: 3, x: 3)>
# array([[ 0.,  4.,  6.],
#        [ 8.,  8., 11.],
#        [13., 15., 15.]], dtype=float32)
# Coordinates:
#   * x        (x) float64 4.5 5.5 6.5
#   * y        (y) float64 1.5 2.5 3.5
# Attributes:
#     long_name:     z
#     actual_range:  [ 0. 15.]

Ideally, the newgrid.dims should be ('height', 'width') instead of ('y', 'x').

Are you willing to help implement and maintain this feature? Yes, but happy to let someone take this up too!