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

**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

```python
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)
#
# 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`](https://www.pygmt.org/v0.6.0/api/generated/pygmt.grdhisteq.equalize_grid.html)

```python
newgrid = pygmt.grdhisteq.equalize_grid(grid=grid)
print(newgrid)
#
# 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!