Merge pull request OSGeo#11982 from dbaston/vrt-common-res

BuildVRT: Add -resolution=compatible

Author: rouault

apps/gdalalg_raster_calc.cpp

@@ -183,17 +183,32 @@ UpdateSourceProperties(SourceProperties &out, const std::string &dsn,
         return std::nullopt;
     }
 
-    // Choose the finest resolution
-    // TODO modify to use common resolution (https://github.com/OSGeo/gdal/issues/11497)
+    // Find a common resolution
     if (source.nX > out.nX)
     {
-        out.nX = source.nX;
-        out.gt[1] = source.gt[1];
+        auto dx = CPLGreatestCommonDivisor(out.gt[1], source.gt[1]);
+        if (dx == 0)
+        {
+            CPLError(CE_Failure, CPLE_AppDefined,
+                     "Failed to find common resolution for inputs.");
+            return std::nullopt;
+        }
+        out.nX = static_cast<int>(
+            std::round(static_cast<double>(out.nX) * out.gt[1] / dx));
+        out.gt[1] = dx;
     }
     if (source.nY > out.nY)
     {
-        out.nY = source.nY;
-        out.gt[5] = source.gt[5];
+        auto dy = CPLGreatestCommonDivisor(out.gt[5], source.gt[5]);
+        if (dy == 0)
+        {
+            CPLError(CE_Failure, CPLE_AppDefined,
+                     "Failed to find common resolution for inputs.");
+            return std::nullopt;
+        }
+        out.nY = static_cast<int>(
+            std::round(static_cast<double>(out.nY) * out.gt[5] / dy));
+        out.gt[5] = dy;
     }
 
     if (srsMismatch)

apps/gdalalg_raster_mosaic.cpp

@@ -54,13 +54,13 @@ GDALRasterMosaicAlgorithm::GDALRasterMosaicAlgorithm()
                    _("Target resolution (in destination CRS units)"),
                    &m_resolution)
                 .SetDefault("same")
-                .SetMetaVar("<xres>,<yres>|same|average|highest|lowest");
+                .SetMetaVar("<xres>,<yres>|same|average|common|highest|lowest");
         arg.AddValidationAction(
             [this, &arg]()
             {
                 const std::string val = arg.Get<std::string>();
                 if (val != "average" && val != "highest" && val != "lowest" &&
-                    val != "same")
+                    val != "same" && val != "common")
                 {
                     const auto aosTokens =
                         CPLStringList(CSLTokenizeString2(val.c_str(), ",", 0));
@@ -70,10 +70,11 @@ GDALRasterMosaicAlgorithm::GDALRasterMosaicAlgorithm()
                         CPLAtof(aosTokens[0]) <= 0 ||
                         CPLAtof(aosTokens[1]) <= 0)
                     {
-                        ReportError(CE_Failure, CPLE_AppDefined,
-                                    "resolution: two comma separated positive "
-                                    "values should be provided, or 'same', "
-                                    "'average', 'highest' or 'lowest'");
+                        ReportError(
+                            CE_Failure, CPLE_AppDefined,
+                            "resolution: two comma separated positive "
+                            "values should be provided, or 'same', "
+                            "'average', 'common', 'highest' or 'lowest'");
                         return false;
                     }
                 }

apps/gdalalg_raster_stack.cpp

@@ -54,13 +54,13 @@ GDALRasterStackAlgorithm::GDALRasterStackAlgorithm()
                    _("Target resolution (in destination CRS units)"),
                    &m_resolution)
                 .SetDefault("same")
-                .SetMetaVar("<xres>,<yres>|same|average|highest|lowest");
+                .SetMetaVar("<xres>,<yres>|same|average|common|highest|lowest");
         arg.AddValidationAction(
             [this, &arg]()
             {
                 const std::string val = arg.Get<std::string>();
                 if (val != "average" && val != "highest" && val != "lowest" &&
-                    val != "same")
+                    val != "same" && val != "common")
                 {
                     const auto aosTokens =
                         CPLStringList(CSLTokenizeString2(val.c_str(), ",", 0));
@@ -70,10 +70,11 @@ GDALRasterStackAlgorithm::GDALRasterStackAlgorithm()
                         CPLAtof(aosTokens[0]) <= 0 ||
                         CPLAtof(aosTokens[1]) <= 0)
                     {
-                        ReportError(CE_Failure, CPLE_AppDefined,
-                                    "resolution: two comma separated positive "
-                                    "values should be provided, or 'same', "
-                                    "'average', 'highest' or 'lowest'");
+                        ReportError(
+                            CE_Failure, CPLE_AppDefined,
+                            "resolution: two comma-separated positive "
+                            "values should be provided, or 'same', "
+                            "'average', 'common', 'highest' or 'lowest'");
                         return false;
                     }
                 }

apps/gdalbuildvrt_lib.cpp

@@ -34,6 +34,7 @@
 #include "commonutils.h"
 #include "cpl_conv.h"
 #include "cpl_error.h"
+#include "cpl_float.h"
 #include "cpl_progress.h"
 #include "cpl_string.h"
 #include "cpl_vsi.h"
@@ -64,7 +65,8 @@ typedef enum
     HIGHEST_RESOLUTION,
     AVERAGE_RESOLUTION,
     SAME_RESOLUTION,
-    USER_RESOLUTION
+    USER_RESOLUTION,
+    COMMON_RESOLUTION,
 } ResolutionStrategy;
 
 struct DatasetProperty
@@ -1016,6 +1018,21 @@ std::string VRTBuilder::AnalyseRaster(GDALDatasetH hDS,
             // ns_res is negative, the highest resolution is the max value.
             ns_res = std::max(ns_res, padfGeoTransform[GEOTRSFRM_NS_RES]);
         }
+        else if (resolutionStrategy == COMMON_RESOLUTION)
+        {
+            we_res = CPLGreatestCommonDivisor(
+                we_res, padfGeoTransform[GEOTRSFRM_WE_RES]);
+            if (we_res == 0)
+            {
+                return "Failed to get common resolution";
+            }
+            ns_res = CPLGreatestCommonDivisor(
+                ns_res, padfGeoTransform[GEOTRSFRM_NS_RES]);
+            if (ns_res == 0)
+            {
+                return "Failed to get common resolution";
+            }
+        }
         else
         {
             we_res = std::max(we_res, padfGeoTransform[GEOTRSFRM_WE_RES]);
@@ -1942,6 +1959,8 @@ GDALDatasetH GDALBuildVRT(const char *pszDest, int nSrcCount,
         eStrategy = LOWEST_RESOLUTION;
     else if (EQUAL(sOptions.osResolution.c_str(), "same"))
         eStrategy = SAME_RESOLUTION;
+    else if (EQUAL(sOptions.osResolution.c_str(), "common"))
+        eStrategy = COMMON_RESOLUTION;
 
     /* If -srcnodata is specified, use it as the -vrtnodata if the latter is not
      */
@@ -2075,7 +2094,7 @@ GDALBuildVRTOptionsGetParser(GDALBuildVRTOptions *psOptions,
                 "the default value which is 'location'."));
 
     argParser->add_argument("-resolution")
-        .metavar("user|average|highest|lowest|same")
+        .metavar("user|average|common|highest|lowest|same")
         .action(
             [psOptions](const std::string &s)
             {
@@ -2084,7 +2103,8 @@ GDALBuildVRTOptionsGetParser(GDALBuildVRTOptions *psOptions,
                     !EQUAL(psOptions->osResolution.c_str(), "average") &&
                     !EQUAL(psOptions->osResolution.c_str(), "highest") &&
                     !EQUAL(psOptions->osResolution.c_str(), "lowest") &&
-                    !EQUAL(psOptions->osResolution.c_str(), "same"))
+                    !EQUAL(psOptions->osResolution.c_str(), "same") &&
+                    !EQUAL(psOptions->osResolution.c_str(), "common"))
                 {
                     throw std::invalid_argument(
                         CPLSPrintf("Illegal resolution value (%s).",

autotest/cpp/test_cpl.cpp

@@ -21,6 +21,7 @@
 
 #include "cpl_compressor.h"
 #include "cpl_error.h"
+#include "cpl_float.h"
 #include "cpl_hash_set.h"
 #include "cpl_list.h"
 #include "cpl_mask.h"
@@ -5763,4 +5764,36 @@ TEST_F(test_cpl, VSIGlob)
     VSIUnlink(osFilenameRadix.c_str());
 }
 
+TEST_F(test_cpl, CPLGreatestCommonDivisor)
+{
+    CPLErrorStateBackuper state(CPLQuietErrorHandler);
+
+    // These tests serve to document the current behavior.
+    // In some cases the results are dependent on various
+    // hardcoded epsilons and it may be appropriate to change
+    // the expected results.
+
+    EXPECT_EQ(CPLGreatestCommonDivisor(3.0, 5.0), 1.0);
+    EXPECT_EQ(CPLGreatestCommonDivisor(3.0 / 3600, 5.0 / 3600), 1.0 / 3600);
+    EXPECT_EQ(CPLGreatestCommonDivisor(5.0 / 3600, 2.5 / 3600), 2.5 / 3600);
+    EXPECT_EQ(CPLGreatestCommonDivisor(1.0 / 10, 1.0), 1.0 / 10);
+    EXPECT_EQ(CPLGreatestCommonDivisor(1.0 / 17, 1.0 / 13), 1.0 / 221);
+    EXPECT_EQ(CPLGreatestCommonDivisor(1.0 / 17, 1.0 / 3600), 1.0 / 61200);
+
+    // GLO-90 resolutoins
+    EXPECT_EQ(CPLGreatestCommonDivisor(3.0 / 3600, 4.5 / 3600), 1.5 / 3600);
+
+    // WorldDEM resolutions
+    EXPECT_EQ(CPLGreatestCommonDivisor(0.4 / 3600, 0.6 / 3600), 0.2 / 3600);
+    EXPECT_EQ(CPLGreatestCommonDivisor(0.6 / 3600, 0.8 / 3600), 0.2 / 3600);
+
+    EXPECT_EQ(CPLGreatestCommonDivisor(M_PI, M_PI / 6), M_PI / 6);
+    EXPECT_EQ(CPLGreatestCommonDivisor(M_PI / 5, M_PI / 6),
+              0);  // Ideally we would get M_PI / 30
+
+    EXPECT_EQ(CPLGreatestCommonDivisor(2.999999, 3.0), 0);
+    EXPECT_EQ(CPLGreatestCommonDivisor(2.9999999, 3.0), 0);
+    EXPECT_EQ(CPLGreatestCommonDivisor(2.99999999, 3.0), 2.99999999);
+}
+
 }  // namespace

autotest/utilities/test_gdalalg_raster_calc.py

@@ -11,6 +11,8 @@
 # SPDX-License-Identifier: MIT
 ###############################################################################
 
+import functools
+import math
 import re
 
 import gdaltest
@@ -281,7 +283,7 @@ def test_gdalalg_calc_different_resolutions(calc, tmp_vsimem):
 
     xmax = 60
     ymax = 60
-    resolutions = [10, 20, 60]
+    resolutions = [30, 20, 60]
 
     inputs = [tmp_vsimem / f"in_{i}.tif" for i in range(len(resolutions))]
     outfile = tmp_vsimem / "out.tif"
@@ -305,8 +307,8 @@ def test_gdalalg_calc_different_resolutions(calc, tmp_vsimem):
     assert calc.Run()
 
     with gdal.Open(outfile) as ds:
-        assert ds.RasterXSize == xmax / min(resolutions)
-        assert ds.RasterYSize == ymax / min(resolutions)
+        assert ds.RasterXSize == xmax / functools.reduce(math.gcd, resolutions)
+        assert ds.RasterYSize == ymax / functools.reduce(math.gcd, resolutions)
 
         assert np.all(ds.ReadAsArray() == sum(resolutions))
 

autotest/utilities/test_gdalalg_raster_mosaic.py

@@ -156,6 +156,28 @@ def test_gdalalg_raster_mosaic_resolution_lowest():
     assert ds.GetGeoTransform() == pytest.approx((2.0, 1.0, 0.0, 49.0, 0.0, -1.0))
 
 
+def test_gdalalg_raster_mosaic_resolution_common():
+
+    # resolution = 3
+    src1_ds = gdal.GetDriverByName("MEM").Create("", 5, 5)
+    src1_ds.SetGeoTransform([2, 3, 0, 49, 0, -3])
+
+    # resolution = 5
+    src2_ds = gdal.GetDriverByName("MEM").Create("", 3, 3)
+    src2_ds.SetGeoTransform([17, 5, 0, 49, 0, -5])
+
+    alg = get_mosaic_alg()
+    alg["input"] = [src1_ds, src2_ds]
+    alg["output"] = ""
+    alg["resolution"] = "common"
+    assert alg.Run()
+    ds = alg["output"].GetDataset()
+
+    assert ds.RasterXSize == 30
+    assert ds.RasterYSize == 15
+    assert ds.GetGeoTransform() == pytest.approx((2.0, 1.0, 0.0, 49.0, 0.0, -1.0))
+
+
 def test_gdalalg_raster_mosaic_resolution_custom():
 
     src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
@@ -215,19 +237,19 @@ def test_gdalalg_raster_mosaic_resolution_invalid():
     alg = get_mosaic_alg()
     with pytest.raises(
         Exception,
-        match="resolution: two comma separated positive values should be provided, or 'same', 'average', 'highest' or 'lowest'",
+        match="resolution: two comma separated positive values should be provided, or ",
     ):
         alg["resolution"] = "invalid"
 
     with pytest.raises(
         Exception,
-        match="resolution: two comma separated positive values should be provided, or 'same', 'average', 'highest' or 'lowest'",
+        match="resolution: two comma separated positive values should be provided, or ",
     ):
         alg["resolution"] = "0.5"
 
     with pytest.raises(
         Exception,
-        match="resolution: two comma separated positive values should be provided, or 'same', 'average', 'highest' or 'lowest'",
+        match="resolution: two comma separated positive values should be provided, or ",
     ):
         alg["resolution"] = "-0.5,-0.5"
 

autotest/utilities/test_gdalalg_raster_stack.py

@@ -0,0 +1,43 @@
+#!/usr/bin/env pytest
+# -*- coding: utf-8 -*-
+###############################################################################
+# Project:  GDAL/OGR Test Suite
+# Purpose:  'gdal raster stack' testing
+# Author:  Daniel Baston
+#
+###############################################################################
+# Copyright (c) 2025, ISciences LLC
+#
+# SPDX-License-Identifier: MIT
+###############################################################################
+
+import pytest
+
+from osgeo import gdal
+
+
+@pytest.fixture()
+def stack():
+    return gdal.GetGlobalAlgorithmRegistry()["raster"]["stack"]
+
+
+def test_gdalalg_raster_stack_resolution_common(stack):
+
+    # resolution = 3
+    src1_ds = gdal.GetDriverByName("MEM").Create("", 5, 5)
+    src1_ds.SetGeoTransform([2, 3, 0, 49, 0, -3])
+
+    # resolution = 5
+    src2_ds = gdal.GetDriverByName("MEM").Create("", 3, 3)
+    src2_ds.SetGeoTransform([2, 5, 0, 49, 0, -5])
+
+    stack["input"] = [src1_ds, src2_ds]
+    stack["output"] = ""
+    stack["resolution"] = "common"
+    assert stack.Run()
+    ds = stack["output"].GetDataset()
+
+    assert ds.RasterCount == 2
+    assert ds.RasterXSize == 15
+    assert ds.RasterYSize == 15
+    assert ds.GetGeoTransform() == pytest.approx((2.0, 1.0, 0.0, 49.0, 0.0, -1.0))

autotest/utilities/test_gdalbuildvrt_lib.py

@@ -923,3 +923,49 @@ def test_gdalbuildvrt_lib_nodata_invalid(tmp_vsimem, dtype, nodata):
         gdal.CE_Warning, "cannot represent the specified NoData value"
     ):
         gdal.BuildVRT("", [tmp_vsimem / "in.tif"], VRTNodata=nodata)
+
+
+###############################################################################
+
+
+# Test --resolution=common
+# Also tested by C++ test_cpl.CPLGreatestCommonDivisor
+@gdaltest.enable_exceptions()
+@pytest.mark.parametrize(
+    "resolutions,expected",
+    [
+        ([5 / 3600, 3 / 3600], 1 / 3600),
+        ([5, 3], 1),
+        ([5 / 3600, 2.5 / 3600], 2.5 / 3600),
+        ([1 / 10, 1], 1 / 10),
+        ([1 / 10, 1 / 3], 1 / 30),
+        ([1 / 17, 1 / 13], 1 / 221),
+        ([1 / 17, 1 / 3600], 1 / 61200),
+        ([2.9999999, 3], "common resolution"),
+    ],
+)
+def test_gdalbuildvrt_resolution_common(tmp_vsimem, resolutions, expected):
+
+    inputs = []
+
+    width = 5
+    height = 5
+
+    for i, res in enumerate(resolutions):
+        fname = tmp_vsimem / f"in_{i}.tif"
+        inputs.append(fname)
+
+        nx = round(width / res)
+        ny = round(height / res)
+
+        with gdal.GetDriverByName("GTiff").Create(fname, nx, ny) as ds:
+            ds.SetGeoTransform((0, res, 0, height, 0, -res))
+
+    if type(expected) is str:
+        with pytest.raises(Exception, match=expected):
+            gdal.BuildVRT("", inputs, resolution="common", strict=True)
+    else:
+        with gdal.BuildVRT("", inputs, resolution="common", strict=True) as ds:
+            gt = ds.GetGeoTransform()
+            assert gt[1] == expected
+            assert -gt[5] == expected