Add transform_precision parameter to merge()

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,

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