reproject: support 3-D rasters on dask, cupy, and dask+cupy (#2027)

[brendancol]

Closes #2027.

Summary

reproject() accepted 3-D (y, x, band) inputs via _validate_raster(ndim=(2, 3)), but only the eager numpy path actually handled them. The dask path advertised a 2-D template for chunks that returned 3-D arrays (lazy shape disagreed with .compute() output), the cupy path tripped a CUDA signature mismatch, and merge() crashed at output DataArray construction.

• _reproject_chunk_cupy: per-band loop matching _reproject_chunk_numpy
• _reproject_dask / _reproject_block_adapter: 3-D template with (out_h, out_w, n_bands) shape so lazy metadata matches actual chunk output
• _reproject_dask_cupy: route 3-D sources through map_blocks (fast path is 2-D-only)
• Empty-chunk early returns in both chunk functions now carry the band axis
• merge(): tighten validator to ndim=(2,) so callers see a clean error instead of a ValueError: different number of dimensions on data and dims crash

Backend coverage

• numpy (was already working; kept under regression tests)
• dask + numpy
• cupy
• dask + cupy

Test plan

• TestReproject3DBackends: lazy shape, computed shape, eager-vs-dask parity, integer dtype round-trip, cupy parity, dask+cupy compute
• TestMerge3DRejection: clean ValueError instead of dim mismatch
• Full test_reproject.py suite still green (267 tests)

[brendancol]

PR Review (self-review after follow-ups in 2b8c725)

Blockers

• None.

Suggestions (addressed in 2b8c725)

• _reproject_chunk_cupy 3-D branch was running cp.where(band_data == nodata, cp.nan, band_data) even on the _use_native_cuda=True path. The 2-D path skips that conversion because the native CUDA kernels handle the sentinel internally. The 3-D branch now mirrors that contract: pre-convert only on the CPU-coords fallback.
• 3-D cupy coverage with nodata=0 (uint8 sentinel) was missing — only the nodata=np.nan path was tested. Added test_reproject_3d_cupy_uint8_sentinel_nodata.

Nits (addressed in 2b8c725)

• Dropped getattr(source_data, 'ndim', 2) defensives in both chunk functions; numpy/cupy/dask all expose ndim directly.
• Added a comment on the 3-D dask template explaining that the band axis is intentionally collapsed to one output chunk regardless of the source's band chunking.

What looks good

• Output rank now matches input rank on every backend; lazy result.shape agrees with result.compute().shape.
• Empty-chunk early returns in both chunk functions carry the band axis, so partial-overlap projections on 3-D sources no longer trip the map_blocks template assembly.
• _reproject_dask_cupy delegates 3-D to map_blocks instead of retrofitting the inline fast path — keeps the diff focused.
• merge() validator tightened to ndim=(2,) — _place_same_crs, _merge_arrays_numpy, and the final xr.DataArray(..., dims=[ydim, xdim]) are all 2-D-only, so the previous ndim=(2, 3) was a promise the implementation never kept.
• Tests cover numpy / dask / cupy / dask+cupy with eager-vs-dask assert_allclose parity, integer dtype round-trips on CPU and GPU, and the issue's exact reproducers.
• 268 tests pass; no regressions in the existing 2-D paths.

Checklist

• Algorithm matches the existing 2-D per-band loop pattern
• All implemented backends produce consistent results (assert_allclose vs numpy)
• NaN handling is correct
• Edge cases covered: uint8 round-trip on CPU and GPU, no-overlap empty chunks
• Dask chunk boundaries handled correctly
• No premature materialization introduced
• Benchmark — 3-D bench would be nice-to-have but the existing bench_reproject_vs_rioxarray.py is 2-D; not blocking for a correctness fix
• README feature matrix — reproject already listed; backend coverage went from advertised-but-broken to actually-supported for 3-D
• Docstrings updated on _reproject_dask and _reproject_block_adapter