geoid_height_raster drops input attrs and mishandles 3D rasters

[brendancol]

geoid_height_raster() in xrspatial/reproject/_vertical.py has two metadata propagation bugs.

First, it drops input attrs. The output DataArray is built with only {'units': 'metres', 'model': model}, so user-provided crs, res, transform, _FillValue, long_name, scale_factor, etc. are silently dropped. The output grid is identical to the input grid, so georeferencing metadata should pass through unchanged.

Second, it mishandles 3D inputs. For a (y, x, band)-shaped input, the function uses raster.dims[-2:] (i.e. ('x', 'band')) as the output dims and looks up raster.coords['x'] and raster.coords['band'] as the lat/lon arrays. This produces garbage: the geoid is interpolated at the wrong coordinates and the output is labelled with the wrong dims and the wrong shape. The validator accepts 3D inputs (per #1431), so this path is reachable from the public API.

Reproduction

import numpy as np
import xarray as xr
from xrspatial.reproject import geoid_height_raster

# Drop attrs
r = xr.DataArray(
    np.zeros((4, 4)), dims=['y', 'x'],
    coords={'y': [3.0, 2, 1, 0], 'x': [0.0, 1, 2, 3]},
    attrs={'crs': 'EPSG:4326', 'transform': (1, 0, -0.5, 0, -1, 3.5)},
)
print(geoid_height_raster(r).attrs)
# {'units': 'metres', 'model': 'EGM96'} -- crs and transform lost

# 3D garbage
r3 = xr.DataArray(
    np.zeros((4, 4, 3)), dims=['y', 'x', 'band'],
    coords={'y': [3.0, 2, 1, 0], 'x': [0.0, 1, 2, 3], 'band': [1, 2, 3]},
    attrs={'crs': 'EPSG:4326'},
)
out = geoid_height_raster(r3)
print(out.shape, out.dims)
# (4, 3) ('x', 'band') -- should be (4, 4) ('y', 'x') with band stripped

Expected behaviour

• Output attrs carries forward input attrs (crs, res, transform, nodata, _FillValue) and adds units and model.
• For 3D inputs, output is 2D with the spatial dims of the input (band axis dropped) since the geoid undulation is purely a function of position.
• For 2D inputs, behaviour is unchanged.

Suggested fix

In _vertical.py::geoid_height_raster:

• Use _find_spatial_dims to identify the y/x dims regardless of the band layout.
• Carry input attrs forward; layer units and model on top.
• Document the 3D-to-2D reduction in the docstring.