Merge pull request OSGeo#12004 from dbaston/gdal-raster-calc-extent-tol

gdal raster calc: use tolerance for extent check

Author: rouault

apps/gdalalg_raster_calc.cpp

@@ -162,8 +162,11 @@ UpdateSourceProperties(SourceProperties &out, const std::string &dsn,
         double ymin =
             source.gt[3] + source.nX * source.gt[4] + source.nY * source.gt[5];
 
-        // TODO use a tolerance here
-        if (xmax != xmaxOut || ymin != yminOut)
+        // Max allowable extent misalignment, expressed as fraction of a pixel
+        constexpr double EXTENT_RTOL = 1e-3;
+
+        if (std::abs(xmax - xmaxOut) > EXTENT_RTOL * std::abs(source.gt[1]) ||
+            std::abs(ymin - yminOut) > EXTENT_RTOL * std::abs(source.gt[5]))
         {
             extentMismatch = true;
         }

autotest/utilities/test_gdalalg_raster_calc.py

@@ -338,6 +338,46 @@ def test_gdalalg_raster_calc_error_extent_mismatch(calc, tmp_vsimem):
         assert src.GetGeoTransform() == dst.GetGeoTransform()
 
 
+def test_gdalalg_raster_calc_error_extent_within_tolerance(calc, tmp_vsimem):
+
+    input_1 = tmp_vsimem / "in1.tif"
+    input_2 = tmp_vsimem / "in2.tif"
+    outfile = tmp_vsimem / "out.tif"
+
+    with gdal.GetDriverByName("GTIff").Create(input_2, 3600, 3600) as ds:
+        ds.SetGeoTransform(
+            (
+                2.4999861111111112e01,
+                2.7777777777777778e-04,
+                0.0000000000000000e00,
+                8.0000138888888884e01,
+                0.0000000000000000e00,
+                -2.7777777777777778e-04,
+            )
+        )
+    with gdal.GetDriverByName("GTiff").Create(input_1, 3600, 3600) as ds:
+        # this geotransform represents error introduced by writing a dataset with the above
+        # geotransform to netCDF
+        ds.SetGeoTransform(
+            (
+                2.4999861111111112e01,
+                2.7777777777777778e-04,
+                0.0000000000000000e00,
+                8.0000138888888884e01,
+                0.0000000000000000e00,
+                -2.7777777777778173e-04,
+            )
+        )
+
+    calc["input"] = [f"A={input_1}", f"B={input_2}"]
+    calc["output"] = outfile
+    calc["calc"] = ["A+B"]
+    assert calc.Run()
+
+    with gdal.Open(input_1) as src, gdal.Open(outfile) as dst:
+        assert src.GetGeoTransform() == dst.GetGeoTransform()
+
+
 def test_gdalalg_raster_calc_error_crs_mismatch(calc, tmp_vsimem):
 
     input_1 = tmp_vsimem / "in1.tif"