From 7989e377e682cc67c84c2eab073fb708aaa39763 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 4 May 2026 08:08:56 -0700 Subject: [PATCH] Add transform_precision parameter to merge() (#1452) reproject() exposes transform_precision but merge() did not, so users had no way to pick the control-grid subdivision count or request exact per-pixel transforms when merging. The parameter is now plumbed through _merge_inmemory and _merge_dask instead of being hardcoded to 16. max_memory was intentionally not added because merge() has no streaming path; the docstring notes the dask path for very large outputs. Closes #1452 --- xrspatial/reproject/__init__.py | 25 +++++++++++----- xrspatial/tests/test_reproject.py | 50 +++++++++++++++++++++++++++++++ 2 files changed, 68 insertions(+), 7 deletions(-) diff --git a/xrspatial/reproject/__init__.py b/xrspatial/reproject/__init__.py index 16f3c1ad..9b416714 100644 --- a/xrspatial/reproject/__init__.py +++ b/xrspatial/reproject/__init__.py @@ -1340,6 +1340,7 @@ def merge( nodata=None, strategy='first', chunk_size=None, + transform_precision=16, name=None, ): """Merge multiple rasters into a single mosaic. @@ -1366,12 +1367,22 @@ def merge( Merge strategy: 'first', 'last', 'mean', 'max', 'min'. chunk_size : int or (int, int) or None Chunk size for dask output. + transform_precision : int + Control-grid subdivisions for the coordinate transform (default 16). + Higher values increase accuracy at the cost of more pyproj calls. + Set to 0 for exact per-pixel transforms matching GDAL/rasterio. name : str or None Name for the output DataArray. Returns ------- xr.DataArray + + Notes + ----- + There is no streaming in-memory path; for very large output mosaics, + pass dask-backed inputs (or rely on the automatic promotion to the + dask path) so that each output chunk is computed independently. """ if not rasters: raise ValueError("merge(): rasters list must not be empty") @@ -1385,7 +1396,7 @@ def merge( bounds=bounds, width=None, height=None, - transform_precision=None, + transform_precision=transform_precision, func_name='merge', ) @@ -1473,12 +1484,12 @@ def merge( if any_dask: result_data = _merge_dask( raster_infos, tgt_wkt, out_bounds, out_shape, - resampling, nd, strategy, chunk_size, + resampling, nd, strategy, chunk_size, transform_precision, ) else: result_data = _merge_inmemory( raster_infos, tgt_wkt, out_bounds, out_shape, - resampling, nd, strategy, + resampling, nd, strategy, transform_precision, ) y_coords, x_coords = _make_output_coords(out_bounds, out_shape) @@ -1563,7 +1574,7 @@ def _place_same_crs(src_data, src_bounds, src_shape, y_desc, def _merge_inmemory( raster_infos, tgt_wkt, out_bounds, out_shape, - resampling, nodata, strategy, + resampling, nodata, strategy, transform_precision, ): """In-memory merge using numpy. @@ -1591,7 +1602,7 @@ def _merge_inmemory( info['src_bounds'], info['src_shape'], info['y_desc'], info['src_wkt'], tgt_wkt, out_bounds, out_shape, - resampling, nodata, 16, + resampling, nodata, transform_precision, ) arrays.append(reprojected) return _merge_arrays_numpy(arrays, nodata, strategy) @@ -1650,7 +1661,7 @@ def _merge_block_adapter( def _merge_dask( raster_infos, tgt_wkt, out_bounds, out_shape, - resampling, nodata, strategy, chunk_size, + resampling, nodata, strategy, chunk_size, transform_precision, ): """Dask merge backend using ``map_blocks``.""" import functools @@ -1696,7 +1707,7 @@ def _merge_dask( resampling=resampling, nodata=nodata, strategy=strategy, - precision=16, + precision=transform_precision, src_footprints_tgt=footprints, same_crs_list=same_crs_list, ) diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py index 893b2150..ab4c7326 100644 --- a/xrspatial/tests/test_reproject.py +++ b/xrspatial/tests/test_reproject.py @@ -1678,6 +1678,56 @@ def test_merge_bounds_rejected(self): with pytest.raises(ValueError, match="right"): merge([self._raster()], bounds=(10, 0, 0, 10)) + def test_merge_accepts_transform_precision_zero(self): + """``transform_precision=0`` requests exact per-pixel transforms.""" + from xrspatial.reproject import merge + raster = _gradient_raster(h=8, w=8) + result = merge([raster], transform_precision=0, resolution=1.0) + assert result.shape[0] > 0 + assert result.shape[1] > 0 + + def test_merge_accepts_transform_precision_default(self): + """Default ``transform_precision`` (16) leaves merge() callable.""" + from xrspatial.reproject import merge + raster = _gradient_raster(h=8, w=8) + result = merge([raster], resolution=1.0) + assert result.shape[0] > 0 + assert result.shape[1] > 0 + + def test_merge_rejects_negative_transform_precision(self): + from xrspatial.reproject import merge + with pytest.raises(ValueError, match="transform_precision"): + merge([self._raster()], transform_precision=-1) + + def test_merge_rejects_float_transform_precision(self): + from xrspatial.reproject import merge + with pytest.raises(ValueError, match="transform_precision"): + merge([self._raster()], transform_precision=1.5) + + def test_merge_transform_precision_threaded_to_chunks(self): + """precision=0 (exact) and precision=16 should agree on smooth inputs. + + For inputs where the control-grid approximation is already very + close to the per-pixel transform, the two paths should give the + same merged output to floating-point tolerance. + """ + from xrspatial.reproject import merge + # Two adjacent same-CRS gradients in EPSG:4326 reprojected to + # the same CRS: the control grid is dense enough that precision=16 + # and precision=0 produce identical numbers. + a = _gradient_raster(h=16, w=16, x_range=(-5, 0), y_range=(-5, 5)) + b = _gradient_raster(h=16, w=16, x_range=(0, 5), y_range=(-5, 5)) + out16 = merge([a, b], target_crs='EPSG:4326', + resolution=1.0, transform_precision=16) + out0 = merge([a, b], target_crs='EPSG:4326', + resolution=1.0, transform_precision=0) + assert out16.shape == out0.shape + v16 = out16.values + v0 = out0.values + valid = ~np.isnan(v16) & ~np.isnan(v0) + assert valid.any() + np.testing.assert_allclose(v0[valid], v16[valid], rtol=1e-10) + # ===================================================================== # Issue #1435: NaN/Inf rejection in scalar inputs