reproject: handle (band, y, x) 3-D inputs (#2182)

[brendancol]

Closes #2182

Summary

• The per-chunk worker assumed the band axis was trailing: it sliced source_data[r:, c:] and read window.shape[2] for the band count. A rasterio/rioxarray-style (band, y, x) input sliced the band/y axes instead of y/x and either crashed with a band coord-length mismatch or returned wrong-shape data.
• The entry point now transposes 3-D inputs to the canonical (y, x, band) layout before dispatch and transposes the output back to the input's dim order, so a (band, y, x) raster round-trips as (band, y, x).

Backends covered

numpy, dask+numpy, cupy, dask+cupy. New tests run on all four (cupy/dask+cupy paths gated on GPU availability via skipif).

Test plan

• pytest xrspatial/tests/test_reproject.py::TestReproject3DBandFirst -- 9 passed
• pytest xrspatial/tests/test_reproject.py -- 277 passed, no regressions
• Pixel-equality check that (band, y, x) and (y, x, band) inputs produce the same data after a common transpose

[brendancol]

PR Review: reproject handles (band, y, x) 3-D inputs (#2182)

Blockers

None.

Suggestions

None.

Nits

• xrspatial/reproject/__init__.py:670: the _canonical tuple is built with _band_dim_in even when it is None. The guard _band_dim_in is not None and _input_dims != _canonical saves the call, but a tuple containing None is a little odd to read. Building the tuple inside the guard would make the intent clearer.
• xrspatial/tests/test_reproject.py:3996: test_band_first_dask_matches_numpy calls _make_band_first_raster() three times. Caching the raster in a local would save two RNG passes (cosmetic; the test is fast either way).

What looks good

• The fix is localized to the entry point: transpose to canonical (y, x, band) before dispatch, transpose back after the output DataArray is built. Downstream workers and the per-block dask adapter are unchanged, which keeps the blast radius small.
• The set-equality check on the final transpose (set(_input_dims) == set(result.dims)) is a defensive guard against crashing if the band-dim name ever differs between input and output.
• Tests cover all four backends (numpy, dask, cupy, dask+cupy) with both layout-preservation and pixel-equality assertions.
• xr.DataArray.transpose is metadata-only on dask arrays, so no extra materialization is introduced for the band-first dask path.

Checklist

• Algorithm matches reference/paper -- N/A (pure layout fix)
• All backends produce consistent results -- verified by test_band_first_matches_band_last and test_band_first_dask_matches_numpy
• NaN handling is correct -- unchanged from the (y, x, band) path
• Edge cases covered -- integer dtype roundtrip included
• Dask chunk boundaries handled -- lazy transpose preserves chunking semantics
• No premature materialization or unnecessary copies -- transpose is metadata-only
• Benchmark not needed -- no perf-sensitive code touched
• README feature matrix -- N/A (no new API)
• Docstrings -- existing docstring still accurate