reproject: output bounds estimation skips datum shifts, wrong bounds for NAD27 and high-res rasters

[brendancol]

Describe the bug

reproject() computes the output raster grid in _compute_output_grid() (xrspatial/reproject/_grid.py), which transforms the source boundary into the target CRS through _transform_boundary(). That helper prefers the Numba fast path transform_points() (xrspatial/reproject/_projections.py) over pyproj.

transform_points() documents that it intentionally skips datum-shift wrapping, on the grounds that the metre-level error is sub-pixel for boundary estimation. But the per-pixel data path try_numba_transform() does apply a Helmert datum shift for the same CRS pairs. So when source and target differ by a datum shift, the output bounds are computed without the shift while the pixel data is reprojected with it. The bounds and the data disagree.

The error is not always sub-pixel. EPSG:4267 (NAD27) is treated as interchangeable with WGS84/NAD83 in _is_supported_geographic(), but NAD27 to WGS84/NAD83 differs by tens to over a hundred metres in the continental US. On a high-resolution raster that is many pixels of error in the computed grid.

Reproduce

import numpy as np, pyproj
from xrspatial.reproject._projections import transform_points

src = pyproj.CRS.from_epsg(4267)  # NAD27
tgt = pyproj.CRS.from_epsg(3857)  # Web Mercator
xs, ys = np.array([-100.0]), np.array([40.0])  # CONUS

fast = transform_points(src, tgt, xs, ys)
ref = pyproj.Transformer.from_crs(src, tgt, always_xy=True).transform(xs, ys)
# fast x differs from pyproj by ~45 m because the NAD27 shift is skipped

Expected behavior

Output bounds estimation should account for the datum shift, or otherwise guarantee correct bounds for high-resolution and grid-shift-datum cases, instead of silently skipping it.

Additional context

Cited lines: xrspatial/reproject/_projections.py:2115 (the "No datum-shift wrapping" note) and xrspatial/reproject/_grid.py:261 (where _compute_output_grid uses it). Found during a code-review sweep.