diff --git a/xrspatial/reproject/__init__.py b/xrspatial/reproject/__init__.py index c3327faf..8519dcf1 100644 --- a/xrspatial/reproject/__init__.py +++ b/xrspatial/reproject/__init__.py @@ -1385,6 +1385,7 @@ def merge( nodata=None, strategy='first', chunk_size=None, + transform_precision=16, name=None, ): """Merge multiple rasters into a single mosaic. @@ -1411,6 +1412,10 @@ 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. @@ -1426,6 +1431,12 @@ def merge( ``time`` coord) are carried through to the output. Coords aligned to the spatial dims are dropped because their values do not apply to the merged grid. + + 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") @@ -1439,7 +1450,7 @@ def merge( bounds=bounds, width=None, height=None, - transform_precision=None, + transform_precision=transform_precision, func_name='merge', ) @@ -1532,12 +1543,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) @@ -1640,7 +1651,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. @@ -1675,7 +1686,7 @@ def _merge_inmemory( info['src_bounds'], info['src_shape'], info['y_desc'], info['src_wkt'], tgt_wkt, out_bounds, out_shape, - resampling, r_nd, 16, + resampling, r_nd, transform_precision, ) # Canonicalize this raster's sentinel to NaN before the merge so # rasters with different sentinels merge correctly. @@ -1754,7 +1765,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 @@ -1803,7 +1814,7 @@ def _merge_dask( resampling=resampling, nodata=nodata, strategy=strategy, - precision=16, + precision=transform_precision, src_footprints_tgt=footprints, raster_nodata_list=rnodata_list, same_crs_list=same_crs_list, diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py index 632ff995..d7e4b504 100644 --- a/xrspatial/tests/test_reproject.py +++ b/xrspatial/tests/test_reproject.py @@ -1801,6 +1801,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