gdaladdo: reuse previous resampling method (from GTiff RESAMPLING metadata item) if not specifying -r and overview levels if not specifying them

Author: rouault

apps/gdaladdo.cpp

@@ -635,10 +635,9 @@ MAIN_START(nArgc, papszArgv)
     if (nArgc < 1)
         exit(-nArgc);
 
-    const char *pszResampling = "nearest";
+    std::string osResampling;
     const char *pszFilename = nullptr;
-    int anLevels[1024] = {};
-    int nLevelCount = 0;
+    std::vector<int> anLevels;
     int nResultStatus = 0;
     bool bReadOnly = false;
     bool bClean = false;
@@ -679,7 +678,7 @@ MAIN_START(nArgc, papszArgv)
         else if (EQUAL(papszArgv[iArg], "-r"))
         {
             CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
-            pszResampling = papszArgv[++iArg];
+            osResampling = papszArgv[++iArg];
         }
         else if (EQUAL(papszArgv[iArg], "-ro"))
         {
@@ -752,11 +751,10 @@ MAIN_START(nArgc, papszArgv)
         {
             pszFilename = papszArgv[iArg];
         }
-        else if (atoi(papszArgv[iArg]) > 0 &&
-                 static_cast<size_t>(nLevelCount) < CPL_ARRAYSIZE(anLevels))
+        else if (atoi(papszArgv[iArg]) > 0)
         {
-            anLevels[nLevelCount++] = atoi(papszArgv[iArg]);
-            if (anLevels[nLevelCount - 1] == 1)
+            anLevels.push_back(atoi(papszArgv[iArg]));
+            if (anLevels.back() == 1)
             {
                 printf(
                     "Warning: Overview with subsampling factor of 1 requested. "
@@ -813,24 +811,54 @@ MAIN_START(nArgc, papszArgv)
     if (hDataset == nullptr)
         exit(2);
 
+    if (!bClean && osResampling.empty())
+    {
+        auto poDS = GDALDataset::FromHandle(hDataset);
+        if (poDS->GetRasterCount() > 0)
+        {
+            auto poBand = poDS->GetRasterBand(1);
+            if (poBand->GetOverviewCount() > 0)
+            {
+                const char *pszResampling =
+                    poBand->GetOverview(0)->GetMetadataItem("RESAMPLING");
+                if (pszResampling)
+                {
+                    osResampling = pszResampling;
+                    if (pfnProgress == GDALDummyProgress)
+                        CPLDebug("GDAL",
+                                 "Reusing resampling method %s from existing "
+                                 "overview",
+                                 pszResampling);
+                    else
+                        printf("Info: reusing resampling method %s from "
+                               "existing overview.\n",
+                               pszResampling);
+                }
+            }
+        }
+        if (osResampling.empty())
+            osResampling = "nearest";
+    }
+
     /* -------------------------------------------------------------------- */
     /*      Clean overviews.                                                */
     /* -------------------------------------------------------------------- */
     if (bClean)
     {
-        if (GDALBuildOverviews(hDataset, pszResampling, 0, nullptr, 0, nullptr,
+        if (GDALBuildOverviews(hDataset, "NONE", 0, nullptr, 0, nullptr,
                                pfnProgress, pProgressArg) != CE_None)
         {
-            printf("Cleaning overviews failed.\n");
+            fprintf(stderr, "Cleaning overviews failed.\n");
             nResultStatus = 200;
         }
     }
     else if (bPartialRefreshFromSourceTimestamp)
     {
         if (!PartialRefreshFromSourceTimestamp(
-                GDALDataset::FromHandle(hDataset), pszResampling, nLevelCount,
-                anLevels, nBandCount, panBandList, bMinSizeSpecified, nMinSize,
-                pfnProgress, pProgressArg))
+                GDALDataset::FromHandle(hDataset), osResampling.c_str(),
+                static_cast<int>(anLevels.size()), anLevels.data(), nBandCount,
+                panBandList, bMinSizeSpecified, nMinSize, pfnProgress,
+                pProgressArg))
         {
             nResultStatus = 1;
         }
@@ -839,18 +867,20 @@ MAIN_START(nArgc, papszArgv)
     {
         if (!PartialRefreshFromProjWin(
                 GDALDataset::FromHandle(hDataset), dfULX, dfULY, dfLRX, dfLRY,
-                pszResampling, nLevelCount, anLevels, nBandCount, panBandList,
-                bMinSizeSpecified, nMinSize, pfnProgress, pProgressArg))
+                osResampling.c_str(), static_cast<int>(anLevels.size()),
+                anLevels.data(), nBandCount, panBandList, bMinSizeSpecified,
+                nMinSize, pfnProgress, pProgressArg))
         {
             nResultStatus = 1;
         }
     }
     else if (bPartialRefreshFromSourceExtent)
     {
         if (!PartialRefreshFromSourceExtent(
-                GDALDataset::FromHandle(hDataset), aosSources, pszResampling,
-                nLevelCount, anLevels, nBandCount, panBandList,
-                bMinSizeSpecified, nMinSize, pfnProgress, pProgressArg))
+                GDALDataset::FromHandle(hDataset), aosSources,
+                osResampling.c_str(), static_cast<int>(anLevels.size()),
+                anLevels.data(), nBandCount, panBandList, bMinSizeSpecified,
+                nMinSize, pfnProgress, pProgressArg))
         {
             nResultStatus = 1;
         }
@@ -863,7 +893,32 @@ MAIN_START(nArgc, papszArgv)
         /* --------------------------------------------------------------------
          */
 
-        if (nLevelCount == 0)
+        // If no levels are specified, reuse the potentially existing ones.
+        if (anLevels.empty())
+        {
+            auto poDS = GDALDataset::FromHandle(hDataset);
+            if (poDS->GetRasterCount() > 0)
+            {
+                auto poBand = poDS->GetRasterBand(1);
+                const int nExistingCount = poBand->GetOverviewCount();
+                if (nExistingCount > 0)
+                {
+                    for (int iOvr = 0; iOvr < nExistingCount; ++iOvr)
+                    {
+                        auto poOverview = poBand->GetOverview(iOvr);
+                        if (poOverview)
+                        {
+                            const int nOvFactor = GDALComputeOvFactor(
+                                poOverview->GetXSize(), poBand->GetXSize(),
+                                poOverview->GetYSize(), poBand->GetYSize());
+                            anLevels.push_back(nOvFactor);
+                        }
+                    }
+                }
+            }
+        }
+
+        if (anLevels.empty())
         {
             const int nXSize = GDALGetRasterXSize(hDataset);
             const int nYSize = GDALGetRasterYSize(hDataset);
@@ -872,20 +927,21 @@ MAIN_START(nArgc, papszArgv)
                    DIV_ROUND_UP(nYSize, nOvrFactor) > nMinSize)
             {
                 nOvrFactor *= 2;
-                anLevels[nLevelCount++] = nOvrFactor;
+                anLevels.push_back(nOvrFactor);
             }
         }
 
         // Only HFA supports selected layers
         if (nBandCount > 0)
             CPLSetConfigOption("USE_RRD", "YES");
 
-        if (nLevelCount > 0 &&
-            GDALBuildOverviews(hDataset, pszResampling, nLevelCount, anLevels,
-                               nBandCount, panBandList, pfnProgress,
-                               pProgressArg) != CE_None)
+        if (!anLevels.empty() &&
+            GDALBuildOverviews(hDataset, osResampling.c_str(),
+                               static_cast<int>(anLevels.size()),
+                               anLevels.data(), nBandCount, panBandList,
+                               pfnProgress, pProgressArg) != CE_None)
         {
-            printf("Overview building failed.\n");
+            fprintf(stderr, "Overview building failed.\n");
             nResultStatus = 100;
         }
     }

autotest/utilities/test_gdaladdo.py

@@ -324,3 +324,54 @@ def test_gdaladdo_partial_refresh_from_source_extent(gdaladdo_path, tmp_path):
             ovr_data_refreshed[idx] = ovr_data_ori[idx]
     assert ovr_data_refreshed == ovr_data_ori
     ds = None
+
+
+###############################################################################
+# Test reuse of previous resampling method and overview levels
+
+
+@pytest.mark.parametrize("read_only", [True, False])
+def test_gdaladdo_reuse_previous_resampling_and_levels(
+    gdaladdo_path, tmp_path, read_only
+):
+
+    tmpfilename = str(tmp_path / "test.tif")
+
+    gdal.Translate(tmpfilename, "../gcore/data/byte.tif")
+
+    out, err = gdaltest.runexternal_out_and_err(
+        f"{gdaladdo_path} -r average {tmpfilename}"
+        + (" -ro" if read_only else "")
+        + " 2 4"
+    )
+    assert "ERROR" not in err, (out, err)
+
+    ds = gdal.Open(tmpfilename)
+    assert ds.GetRasterBand(1).GetOverview(0).GetMetadataItem("RESAMPLING") == "AVERAGE"
+    assert ds.GetRasterBand(1).GetOverview(0).Checksum() == 1152
+    ds = None
+
+    # Change resampling method to CUBIC
+    out, err = gdaltest.runexternal_out_and_err(
+        f"{gdaladdo_path} -r cubic {tmpfilename}" + (" -ro" if read_only else "")
+    )
+    assert "ERROR" not in err, (out, err)
+
+    ds = gdal.Open(tmpfilename, gdal.GA_Update)
+    assert ds.GetRasterBand(1).GetOverview(0).GetMetadataItem("RESAMPLING") == "CUBIC"
+    assert ds.GetRasterBand(1).GetOverview(0).Checksum() == 1059
+    # Zeroize overview
+    ds.GetRasterBand(1).GetOverview(0).Fill(0)
+    ds = None
+
+    # Invoke gdaladdo without arguments and check overviews are regenerated
+    # using CUBIC
+    out, err = gdaltest.runexternal_out_and_err(
+        f"{gdaladdo_path} {tmpfilename}" + (" -ro" if read_only else "")
+    )
+    assert "ERROR" not in err, (out, err)
+
+    ds = gdal.Open(tmpfilename)
+    assert ds.GetRasterBand(1).GetOverview(0).GetMetadataItem("RESAMPLING") == "CUBIC"
+    assert ds.GetRasterBand(1).GetOverview(0).Checksum() == 1059
+    ds = None

doc/source/programs/gdaladdo.rst

@@ -35,9 +35,16 @@ most supported file formats with one of several downsampling algorithms.
 
 .. option:: -r {nearest|average|rms|gauss|cubic|cubicspline|lanczos|average_magphase|mode}
 
-    Select a resampling algorithm.
+    Select a resampling algorithm. The default is ``nearest``, which is generally not
+    appropriate if sub-pixel accuracy is desired.
 
-    ``nearest`` applies a nearest neighbour (simple sampling) resampler (default)
+    Starting with GDAL 3.9, when refreshing existing TIFF overviews, the previously
+    used method, as noted in the RESAMPLING metadata item of the overview, will
+    be used if :option:`-r` is not specified.
+
+    The available methods are:
+
+    ``nearest`` applies a nearest neighbour (simple sampling) resampler.
 
     ``average`` computes the average of all non-NODATA contributing pixels. Starting with GDAL 3.1, this is a weighted average taking into account properly the weight of source pixels not contributing fully to the target pixel.
 
@@ -131,10 +138,15 @@ most supported file formats with one of several downsampling algorithms.
 
     .. versionadded:: 2.3
 
-        levels are no longer required to build overviews.
+        Levels are no longer required to build overviews.
         In which case, appropriate overview power-of-two factors will be selected
         until the smallest overview is smaller than the value of the -minsize switch.
 
+        Starting with GDAL 3.9, if there are already existing overviews, the
+        corresponding levels will be used to refresh them if no explicit levels
+        are specified.
+
+
 gdaladdo will honour properly NODATA_VALUES tuples (special dataset metadata) so
 that only a given RGB triplet (in case of a RGB image) will be considered as the
 nodata value and not each value of the triplet independently per band.