From f0e05c2a369a146a408de7c27b262419961fa61f Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 18 May 2026 07:37:27 -0700 Subject: [PATCH 1/2] reproject: support 3-D rasters on dask, cupy, and dask+cupy (#2027) The 3-D (y, x, band) input contract was only honoured on the eager numpy path. The dask backend built a 2-D template that disagreed with its own 3-D chunks (lazy DataArray claimed 2-D, .compute() crashed). The cupy backend never gained a 3-D branch and tripped a CUDA signature mismatch. merge() advertised 3-D via the validator but crashed at output DataArray construction. - _reproject_chunk_cupy: add per-band loop matching the numpy path - _reproject_dask: build a 3-D template (out_h, out_w, n_bands) for 3-D sources so the lazy metadata matches actual chunk output - _reproject_dask_cupy: route 3-D sources through map_blocks (the fast path's inline window loop is 2-D-only) - _reproject_block_adapter: read 2-D location from a 3-D template; return matching-rank empty chunks - _reproject_chunk_numpy / _cupy: empty-chunk early returns now carry the band axis when the source is 3-D - merge(): tighten _validate_raster to ndim=(2,) so the API surface stops promising 3-D support the implementation cannot deliver - tests: cross-backend 3-D parity on numpy / dask / cupy / dask+cupy plus a merge 3-D rejection test --- xrspatial/reproject/__init__.py | 132 +++++++++++++++++++--- xrspatial/tests/test_reproject.py | 175 ++++++++++++++++++++++++++++++ 2 files changed, 291 insertions(+), 16 deletions(-) diff --git a/xrspatial/reproject/__init__.py b/xrspatial/reproject/__init__.py index ce2f59dd8..e4d0d4146 100644 --- a/xrspatial/reproject/__init__.py +++ b/xrspatial/reproject/__init__.py @@ -251,10 +251,19 @@ def _reproject_chunk_numpy( c_min = np.nanmin(src_col_px) c_max = np.nanmax(src_col_px) + # 3-D source: empty-chunk returns must carry the band axis or the + # dask map_blocks template (which is 3-D for 3-D sources) sees a + # shape mismatch (#2027). + _src_ndim = getattr(source_data, 'ndim', 2) + if _src_ndim == 3: + _empty_shape = (*chunk_shape, source_data.shape[2]) + else: + _empty_shape = chunk_shape + if not np.isfinite(r_min) or not np.isfinite(r_max): - return np.full(chunk_shape, nodata, dtype=np.float64) + return np.full(_empty_shape, nodata, dtype=np.float64) if not np.isfinite(c_min) or not np.isfinite(c_max): - return np.full(chunk_shape, nodata, dtype=np.float64) + return np.full(_empty_shape, nodata, dtype=np.float64) r_min = int(np.floor(r_min)) - 2 r_max = int(np.ceil(r_max)) + 3 @@ -263,7 +272,7 @@ def _reproject_chunk_numpy( # Check overlap if r_min >= src_h or r_max <= 0 or c_min >= src_w or c_max <= 0: - return np.full(chunk_shape, nodata, dtype=np.float64) + return np.full(_empty_shape, nodata, dtype=np.float64) # Clip to source bounds r_min_clip = max(0, r_min) @@ -276,7 +285,7 @@ def _reproject_chunk_numpy( _MAX_WINDOW_PIXELS = 64 * 1024 * 1024 # 64 Mpix (~512 MB for float64) win_pixels = (r_max_clip - r_min_clip) * (c_max_clip - c_min_clip) if win_pixels > _MAX_WINDOW_PIXELS: - return np.full(chunk_shape, nodata, dtype=np.float64) + return np.full(_empty_shape, nodata, dtype=np.float64) # Extract source window window = source_data[r_min_clip:r_max_clip, c_min_clip:c_max_clip] @@ -337,6 +346,14 @@ def _reproject_chunk_cupy( src_crs = _crs_from_wkt(src_wkt) tgt_crs = _crs_from_wkt(tgt_wkt) + # 3-D source: empty-chunk returns must carry the band axis to match + # the dask+cupy map_blocks template (#2027). + _src_ndim = getattr(source_data, 'ndim', 2) + if _src_ndim == 3: + _empty_shape = (*chunk_shape, source_data.shape[2]) + else: + _empty_shape = chunk_shape + # Try CUDA transform first (keeps coordinates on-device) cuda_result = None if src_crs is not None and tgt_crs is not None: @@ -372,7 +389,7 @@ def _reproject_chunk_cupy( ) if not (np.isfinite(r_min_val) and np.isfinite(r_max_val) and np.isfinite(c_min_val) and np.isfinite(c_max_val)): - return cp.full(chunk_shape, nodata, dtype=cp.float64) + return cp.full(_empty_shape, nodata, dtype=cp.float64) r_min = int(np.floor(r_min_val)) - 2 r_max = int(np.ceil(r_max_val)) + 3 c_min = int(np.floor(c_min_val)) - 2 @@ -407,9 +424,9 @@ def _reproject_chunk_cupy( c_min = np.nanmin(src_col_px) c_max = np.nanmax(src_col_px) if not np.isfinite(r_min) or not np.isfinite(r_max): - return cp.full(chunk_shape, nodata, dtype=cp.float64) + return cp.full(_empty_shape, nodata, dtype=cp.float64) if not np.isfinite(c_min) or not np.isfinite(c_max): - return cp.full(chunk_shape, nodata, dtype=cp.float64) + return cp.full(_empty_shape, nodata, dtype=cp.float64) r_min = int(np.floor(r_min)) - 2 r_max = int(np.ceil(r_max)) + 3 c_min = int(np.floor(c_min)) - 2 @@ -417,7 +434,7 @@ def _reproject_chunk_cupy( _use_native_cuda = False if r_min >= src_h or r_max <= 0 or c_min >= src_w or c_max <= 0: - return cp.full(chunk_shape, nodata, dtype=cp.float64) + return cp.full(_empty_shape, nodata, dtype=cp.float64) r_min_clip = max(0, r_min) r_max_clip = min(src_h, r_max) @@ -429,7 +446,7 @@ def _reproject_chunk_cupy( _MAX_WINDOW_PIXELS = 64 * 1024 * 1024 # 64 Mpix (~512 MB for float64) win_pixels = (r_max_clip - r_min_clip) * (c_max_clip - c_min_clip) if win_pixels > _MAX_WINDOW_PIXELS: - return cp.full(chunk_shape, nodata, dtype=cp.float64) + return cp.full(_empty_shape, nodata, dtype=cp.float64) window = source_data[r_min_clip:r_max_clip, c_min_clip:c_max_clip] if hasattr(window, 'compute'): @@ -437,12 +454,42 @@ def _reproject_chunk_cupy( if not isinstance(window, cp.ndarray): window = cp.asarray(window) orig_dtype = window.dtype - window = window.astype(cp.float64) # Adjust coordinates relative to window (stays on GPU if CuPy) local_row = src_row_px - r_min_clip local_col = src_col_px - c_min_clip + # Multi-band: reproject each band separately, share coordinates. + # Matches the 3-D branch in _reproject_chunk_numpy so 3-D inputs work + # on cupy and dask+cupy backends instead of crashing with a CUDA + # signature mismatch (#2027). + if window.ndim == 3: + n_bands = window.shape[2] + bands = [] + for b in range(n_bands): + band_data = window[:, :, b].astype(cp.float64) + if not np.isnan(nodata): + band_data = cp.where(band_data == nodata, cp.nan, band_data) + if _use_native_cuda: + band_result = _resample_cupy_native( + band_data, local_row, local_col, + resampling=resampling, nodata=nodata, + ) + else: + band_result = _resample_cupy( + band_data, local_row, local_col, + resampling=resampling, nodata=nodata, + ) + if np.issubdtype(orig_dtype, np.integer): + info = np.iinfo(orig_dtype) + band_result = cp.clip( + cp.round(band_result), info.min, info.max, + ).astype(orig_dtype) + bands.append(band_result) + return cp.stack(bands, axis=-1) + + window = window.astype(cp.float64) + if _use_native_cuda: # Coordinates are already CuPy arrays -- use native CUDA kernels # (nodata->NaN conversion is handled inside _resample_cupy_native) @@ -1233,6 +1280,18 @@ def _reproject_dask_cupy( src_res_x = (src_right - src_left) / src_w src_res_y = (src_top - src_bottom) / src_h + # 3-D source: the fast path's inline loop assumes 2-D windows. + # Delegate to the map_blocks path which handles 3-D via + # _reproject_chunk_cupy's per-band loop (#2027). + if raster.data.ndim == 3: + return _reproject_dask( + raster, src_bounds, src_shape, y_desc, + src_wkt, tgt_wkt, + out_bounds, out_shape, + resampling, nodata, precision, + chunk_size or 2048, True, # is_cupy=True + ) + # Memory check: if the full output doesn't fit in GPU memory, # fall back to the map_blocks path which is O(chunk_size) memory. estimated = out_shape[0] * out_shape[1] * 8 # float64 @@ -1444,21 +1503,37 @@ def _reproject_block_adapter( src_wkt, tgt_wkt, out_bounds, out_shape, resampling, nodata, precision, - is_cupy, src_footprint_tgt, + is_cupy, src_footprint_tgt, n_bands=None, ): """``map_blocks`` adapter for reprojection. Derives chunk bounds from *block_info* and delegates to the per-chunk worker. + + For 3-D sources the template carries the band axis, so each block + is ``(rh, rw, n_bands)``. The adapter strips the trailing band axis + when computing 2-D chunk bounds and the per-chunk worker returns a + 3-D result that fits the template (#2027). """ info = block_info[0] - (row_start, row_end), (col_start, col_end) = info['array-location'] + # 3-D template -> array-location is 3 entries; spatial dims are the + # first two. Band dim spans the full axis (single chunk). + spatial_loc = info['array-location'][:2] + (row_start, row_end), (col_start, col_end) = spatial_loc chunk_shape = (row_end - row_start, col_end - col_start) cb = _chunk_bounds(out_bounds, out_shape, row_start, row_end, col_start, col_end) + is_3d = n_bands is not None + # Skip chunks that don't overlap the source footprint if src_footprint_tgt is not None and not _bounds_overlap(cb, src_footprint_tgt): + if is_3d: + empty_shape = (*chunk_shape, n_bands) + if is_cupy: + import cupy as cp + return cp.full(empty_shape, nodata, dtype=cp.float64) + return np.full(empty_shape, nodata, dtype=np.float64) return np.full(chunk_shape, nodata, dtype=np.float64) chunk_fn = _reproject_chunk_cupy if is_cupy else _reproject_chunk_numpy @@ -1487,6 +1562,13 @@ def _reproject_dask( adding the full source as a dependency of every output block (which would cause a MemoryError on distributed schedulers when the source exceeds the worker memory limit). + + For 3-D sources with shape ``(H, W, n_bands)`` the template is built + as ``(out_H, out_W, n_bands)`` so the lazy metadata matches the + actual chunk output shape (#2027). Without this, the lazy DataArray + advertised 2-D shape while the underlying chunks produced 3-D + arrays, causing a ``ValueError: replacement data has shape ...`` + crash on ``.compute()``. """ import functools @@ -1499,6 +1581,11 @@ def _reproject_dask( src_bounds, src_wkt, tgt_wkt ) + # Detect 3-D source: chunks will return (rh, rw, n_bands) so the + # template must carry the band axis through to the lazy DataArray. + is_3d = raster.data.ndim == 3 + n_bands = raster.data.shape[2] if is_3d else None + # Bind the source dask array and all scalar params via partial so # map_blocks doesn't detect them as dask Array kwargs (which would # add the full source as a dependency of every output block). @@ -1517,6 +1604,7 @@ def _reproject_dask( precision=precision, is_cupy=is_cupy, src_footprint_tgt=src_footprint_tgt, + n_bands=n_bands, ) # Pick the template dtype to match the eager path: integer sources @@ -1529,9 +1617,16 @@ def _reproject_dask( else: out_dtype = np.dtype(np.float64) - template = da.empty( - out_shape, dtype=out_dtype, chunks=(row_chunks, col_chunks) - ) + if is_3d: + template = da.empty( + (*out_shape, n_bands), + dtype=out_dtype, + chunks=(row_chunks, col_chunks, (n_bands,)), + ) + else: + template = da.empty( + out_shape, dtype=out_dtype, chunks=(row_chunks, col_chunks) + ) return da.map_blocks( bound_adapter, @@ -1635,8 +1730,13 @@ def merge( raise ValueError("merge(): rasters list must not be empty") for i, r in enumerate(rasters): + # merge() only supports 2-D rasters: the merge strategies, same-CRS + # placement, and output DataArray construction all assume (y, x). + # The 3-D (y, x, band) path was never implemented end-to-end and + # crashed at DataArray construction with a dims-vs-shape mismatch + # (#2027). Reject 3-D up front so callers get a clear error. _validate_raster(r, func_name='merge', name=f'rasters[{i}]', - ndim=(2, 3)) + ndim=(2,)) _validate_grid_params( resolution=resolution, diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py index 058e1762f..7356d029c 100644 --- a/xrspatial/tests/test_reproject.py +++ b/xrspatial/tests/test_reproject.py @@ -3694,3 +3694,178 @@ def test_reproject_preserves_lat_lon_dim_names_dask(self): raster.data = da.from_array(raster.values, chunks=(4, 4)) result = reproject(raster, 'EPSG:3857', chunk_size=4) assert result.dims == ('lat', 'lon') + + +# ===================================================================== +# Issue #2027: 3-D (y, x, band) inputs across all backends +# ===================================================================== + +class TestReproject3DBackends: + """reproject() must honour the band axis on every backend. + + The 2-D path worked for years; the dask, cupy, and dask+cupy paths + either silently dropped the band dim from the lazy DataArray or + crashed with a CUDA signature mismatch on 3-D inputs (#2027). + """ + + @staticmethod + def _make_3d_raster(rng_seed=0, h=32, w=32, n_bands=3, dtype=np.float32): + rng = np.random.default_rng(rng_seed) + data = rng.random((h, w, n_bands), dtype=np.float32).astype(dtype) + return xr.DataArray( + data, + dims=['y', 'x', 'band'], + coords={ + 'y': np.linspace(55, 45, h), + 'x': np.linspace(-5, 5, w), + 'band': list(range(n_bands)), + }, + attrs={'crs': 'EPSG:4326', 'nodata': np.nan}, + ) + + def test_reproject_3d_numpy(self): + """Baseline: 3-D numpy reproject keeps band dim.""" + from xrspatial.reproject import reproject + raster = self._make_3d_raster() + result = reproject(raster, 'EPSG:32633') + assert result.ndim == 3 + assert result.dims == ('y', 'x', 'band') + assert result.shape[2] == 3 + # Computed values should be finite for at least part of the output + assert np.any(np.isfinite(result.values)) + + @pytest.mark.skipif(not HAS_DASK, reason="dask required") + def test_reproject_3d_dask_lazy_shape(self): + """Lazy dask DataArray must advertise 3-D shape (not 2-D).""" + from xrspatial.reproject import reproject + raster = self._make_3d_raster() + raster = raster.copy( + data=da.from_array(raster.values, chunks=(16, 16, 3)) + ) + result = reproject(raster, 'EPSG:32633') + assert result.ndim == 3 + assert result.dims == ('y', 'x', 'band') + assert result.shape[2] == 3 + + @pytest.mark.skipif(not HAS_DASK, reason="dask required") + def test_reproject_3d_dask_compute(self): + """Computed dask result keeps band axis without ValueError.""" + from xrspatial.reproject import reproject + raster = self._make_3d_raster() + raster = raster.copy( + data=da.from_array(raster.values, chunks=(16, 16, 3)) + ) + result = reproject(raster, 'EPSG:32633').compute() + assert result.ndim == 3 + assert result.shape[2] == 3 + assert np.any(np.isfinite(result.values)) + + @pytest.mark.skipif(not HAS_DASK, reason="dask required") + def test_reproject_3d_dask_matches_numpy(self): + """Dask 3-D output should match eager numpy output pixel-for-pixel.""" + from xrspatial.reproject import reproject + raster = self._make_3d_raster() + eager = reproject(raster, 'EPSG:32633') + lazy_src = raster.copy( + data=da.from_array(raster.values, chunks=(16, 16, 3)) + ) + lazy = reproject(lazy_src, 'EPSG:32633').compute() + np.testing.assert_allclose( + np.asarray(eager.values), np.asarray(lazy.values), + rtol=1e-6, atol=1e-6, equal_nan=True, + ) + + @pytest.mark.skipif(not HAS_DASK, reason="dask required") + def test_reproject_3d_dask_uint8_dtype_roundtrip(self): + """Integer 3-D dask inputs round-trip to source dtype.""" + from xrspatial.reproject import reproject + rng = np.random.default_rng(1) + data = rng.integers(0, 255, (32, 32, 3), dtype=np.uint8) + raster = xr.DataArray( + da.from_array(data, chunks=(16, 16, 3)), + dims=['y', 'x', 'band'], + coords={'y': np.linspace(55, 45, 32), 'x': np.linspace(-5, 5, 32)}, + attrs={'crs': 'EPSG:4326', 'nodata': 0}, + ) + result = reproject(raster, 'EPSG:32633').compute() + assert result.dtype == np.uint8 + assert result.shape[2] == 3 + + @pytest.mark.skipif(not HAS_CUPY, reason="CuPy not installed") + def test_reproject_3d_cupy(self): + """CuPy 3-D reproject keeps band dim without CUDA signature crash.""" + from xrspatial.reproject import reproject + host = self._make_3d_raster() + gpu_data = cp.asarray(host.values) + raster = host.copy(data=gpu_data) + result = reproject(raster, 'EPSG:32633') + assert result.ndim == 3 + assert result.shape[2] == 3 + # Pull back to host to verify finite values + out = cp.asnumpy(result.data) if isinstance(result.data, cp.ndarray) \ + else np.asarray(result.values) + assert np.any(np.isfinite(out)) + + @pytest.mark.skipif(not HAS_CUPY, reason="CuPy not installed") + def test_reproject_3d_cupy_matches_numpy(self): + from xrspatial.reproject import reproject + host = self._make_3d_raster() + eager = reproject(host, 'EPSG:32633').values + gpu = host.copy(data=cp.asarray(host.values)) + gpu_out = reproject(gpu, 'EPSG:32633') + gpu_arr = cp.asnumpy(gpu_out.data) if isinstance(gpu_out.data, cp.ndarray) \ + else np.asarray(gpu_out.values) + np.testing.assert_allclose( + eager, gpu_arr, rtol=1e-4, atol=1e-4, equal_nan=True, + ) + + @pytest.mark.skipif( + not (HAS_CUPY and HAS_DASK), reason="CuPy + dask required", + ) + def test_reproject_3d_dask_cupy(self): + """dask+cupy 3-D reproject keeps band dim.""" + from xrspatial.reproject import reproject + host = self._make_3d_raster() + gpu_data = da.from_array(cp.asarray(host.values), chunks=(16, 16, 3)) + raster = host.copy(data=gpu_data) + result = reproject(raster, 'EPSG:32633') + assert result.ndim == 3 + assert result.dims == ('y', 'x', 'band') + computed = result.compute() + assert computed.shape[2] == 3 + + +@pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not installed") +class TestMerge3DRejection: + """merge() must reject 3-D inputs with a clear error (#2027). + + Before the fix, merge() advertised 3-D support via its validator but + crashed at output DataArray construction because the merge strategies, + same-CRS placement, and final `dims=[ydim, xdim]` all assume 2-D. We + tighten the validator so callers see a clean message instead. + """ + + def test_merge_rejects_3d_dataarray(self): + from xrspatial.reproject import merge + a = xr.DataArray( + np.random.rand(8, 8, 3), + dims=['y', 'x', 'band'], + coords={ + 'y': np.linspace(5, -5, 8), + 'x': np.linspace(-5, 0, 8), + 'band': [1, 2, 3], + }, + attrs={'crs': 'EPSG:4326'}, + ) + b = xr.DataArray( + np.random.rand(8, 8, 3), + dims=['y', 'x', 'band'], + coords={ + 'y': np.linspace(5, -5, 8), + 'x': np.linspace(0, 5, 8), + 'band': [1, 2, 3], + }, + attrs={'crs': 'EPSG:4326'}, + ) + with pytest.raises(ValueError, match=r"must be 2D"): + merge([a, b], resolution=1.0) From 2b8c7255895ed2686407baafcaebc6b47116313f Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 18 May 2026 07:53:39 -0700 Subject: [PATCH 2/2] Address review nits: cupy nodata, defensives, band-chunk doc (#2027) - _reproject_chunk_cupy: skip the band-data nodata->NaN pre-conversion on the native CUDA path so the 3-D branch matches the 2-D path's contract (the CUDA kernels handle the sentinel internally). - Replace defensive getattr(source_data, 'ndim', 2) shims with direct attribute access; numpy/cupy/dask arrays all expose ndim. - Document why the dask template collapses the band axis to one chunk. - Add cupy uint8 sentinel-nodata round-trip test that exercises the non-NaN nodata branch on the GPU. --- xrspatial/reproject/__init__.py | 22 ++++++++++++++++------ xrspatial/tests/test_reproject.py | 19 +++++++++++++++++++ 2 files changed, 35 insertions(+), 6 deletions(-) diff --git a/xrspatial/reproject/__init__.py b/xrspatial/reproject/__init__.py index e4d0d4146..cfab24881 100644 --- a/xrspatial/reproject/__init__.py +++ b/xrspatial/reproject/__init__.py @@ -254,8 +254,7 @@ def _reproject_chunk_numpy( # 3-D source: empty-chunk returns must carry the band axis or the # dask map_blocks template (which is 3-D for 3-D sources) sees a # shape mismatch (#2027). - _src_ndim = getattr(source_data, 'ndim', 2) - if _src_ndim == 3: + if source_data.ndim == 3: _empty_shape = (*chunk_shape, source_data.shape[2]) else: _empty_shape = chunk_shape @@ -348,8 +347,7 @@ def _reproject_chunk_cupy( # 3-D source: empty-chunk returns must carry the band axis to match # the dask+cupy map_blocks template (#2027). - _src_ndim = getattr(source_data, 'ndim', 2) - if _src_ndim == 3: + if source_data.ndim == 3: _empty_shape = (*chunk_shape, source_data.shape[2]) else: _empty_shape = chunk_shape @@ -468,14 +466,20 @@ def _reproject_chunk_cupy( bands = [] for b in range(n_bands): band_data = window[:, :, b].astype(cp.float64) - if not np.isnan(nodata): - band_data = cp.where(band_data == nodata, cp.nan, band_data) if _use_native_cuda: + # Native CUDA kernels do the nodata->NaN conversion + # internally; matching the 2-D path above. band_result = _resample_cupy_native( band_data, local_row, local_col, resampling=resampling, nodata=nodata, ) else: + # CPU coords path needs explicit conversion before + # cupyx.scipy.ndimage.map_coordinates. + if not np.isnan(nodata): + band_data = cp.where( + band_data == nodata, cp.nan, band_data, + ) band_result = _resample_cupy( band_data, local_row, local_col, resampling=resampling, nodata=nodata, @@ -1618,6 +1622,12 @@ def _reproject_dask( out_dtype = np.dtype(np.float64) if is_3d: + # Band axis is one chunk in the template regardless of how the + # source dask array is chunked along its band axis. The per-block + # worker reads all bands together (via a 2-D y/x slice that + # rejoins band chunks on compute) and emits the full band stack + # for its (rh, rw) tile, so multi-chunk output along bands would + # never get filled. template = da.empty( (*out_shape, n_bands), dtype=out_dtype, diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py index 7356d029c..6edd7012f 100644 --- a/xrspatial/tests/test_reproject.py +++ b/xrspatial/tests/test_reproject.py @@ -3834,6 +3834,25 @@ def test_reproject_3d_dask_cupy(self): computed = result.compute() assert computed.shape[2] == 3 + @pytest.mark.skipif(not HAS_CUPY, reason="CuPy not installed") + def test_reproject_3d_cupy_uint8_sentinel_nodata(self): + """3-D cupy with integer sentinel nodata round-trips to source dtype. + + Exercises the non-NaN nodata path that the float tests skip. + """ + from xrspatial.reproject import reproject + rng = np.random.default_rng(2) + host = rng.integers(0, 255, (32, 32, 3), dtype=np.uint8) + raster = xr.DataArray( + cp.asarray(host), + dims=['y', 'x', 'band'], + coords={'y': np.linspace(55, 45, 32), 'x': np.linspace(-5, 5, 32)}, + attrs={'crs': 'EPSG:4326', 'nodata': 0}, + ) + result = reproject(raster, 'EPSG:32633') + assert result.dtype == np.uint8 + assert result.shape[2] == 3 + @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not installed") class TestMerge3DRejection: