Merge pull request #10602 from rouault/fix_10600

gdal_translate/GDALOverviewDataset: fix half-pixel shift issue when rescaling RPC

apps/gdal_translate_lib.cpp

@@ -1826,10 +1826,24 @@ GDALDatasetH GDALTranslate(const char *pszDest, GDALDatasetH hSrcDataset,
 
             const double df2 = adfSrcWinOri[2];
             const double df3 = adfSrcWinOri[3];
-            dfSAMP_OFF *= nOXSize / df2;
-            dfLINE_OFF *= nOYSize / df3;
-            dfSAMP_SCALE *= nOXSize / df2;
-            dfLINE_SCALE *= nOYSize / df3;
+            const double dfXRatio = nOXSize / df2;
+            const double dfYRatio = nOYSize / df3;
+
+            // For line offset and pixel offset, we need to convert from RPC
+            // pixel center registration convention to GDAL pixel top-left corner
+            // registration convention by adding an initial 0.5 shift, and un-apply
+            // it after scaling.
+
+            dfSAMP_OFF += 0.5;
+            dfSAMP_OFF *= dfXRatio;
+            dfSAMP_OFF -= 0.5;
+
+            dfLINE_OFF += 0.5;
+            dfLINE_OFF *= dfYRatio;
+            dfLINE_OFF -= 0.5;
+
+            dfSAMP_SCALE *= dfXRatio;
+            dfLINE_SCALE *= dfYRatio;
 
             CPLString osField;
             osField.Printf("%.15g", dfLINE_OFF);

autotest/gcore/overviewds.py

@@ -225,13 +225,12 @@ def test_overviewds_4(tmp_path):
 
     for key in rpc_md:
         assert ds.GetMetadataItem(key, "RPC") == got_md[key]
-        if (
-            key == "LINE_SCALE"
-            or key == "SAMP_SCALE"
-            or key == "LINE_OFF"
-            or key == "SAMP_OFF"
-        ):
-            assert float(got_md[key]) == myfloat(rpc_md[key]) / 2
+        if key == "LINE_SCALE" or key == "SAMP_SCALE":
+            assert float(got_md[key]) == pytest.approx(myfloat(rpc_md[key]) / 2)
+        elif key == "LINE_OFF" or key == "SAMP_OFF":
+            assert float(got_md[key]) == pytest.approx(
+                (myfloat(rpc_md[key]) + 0.5) / 2 - 0.5
+            )
         elif got_md[key] != rpc_md[key]:
             print(got_md[key])
             print(rpc_md[key])

autotest/utilities/test_gdal_translate.py

@@ -872,16 +872,32 @@ def test_gdal_translate_32(gdal_translate_path, tmp_path):
 
     dst_tif = str(tmp_path / "out.tif")
 
+    src_ds = gdal.Open("../gcore/data/byte_rpc.tif")
+    src_md = src_ds.GetMetadata("RPC")
+    srcxoff = 1
+    srcyoff = 2
+    srcwidth = 13
+    srcheight = 14
+    widthratio = 200
+    heightratio = 300
     gdaltest.runexternal(
-        f"{gdal_translate_path} ../gcore/data/byte_rpc.tif {dst_tif} -srcwin 1 2 13 14 -outsize 150% 300%"
+        f"{gdal_translate_path} ../gcore/data/byte_rpc.tif {dst_tif} -srcwin {srcxoff} {srcyoff} {srcwidth} {srcheight} -outsize {widthratio}% {heightratio}%"
     )
+    widthratio /= 100.0
+    heightratio /= 100.0
     ds = gdal.Open(dst_tif)
     md = ds.GetMetadata("RPC")
-    assert (
-        float(md["LINE_OFF"]) == pytest.approx(47496, abs=1e-5)
-        and float(md["LINE_SCALE"]) == pytest.approx(47502, abs=1e-5)
-        and float(md["SAMP_OFF"]) == pytest.approx(19676.6923076923, abs=1e-5)
-        and float(md["SAMP_SCALE"]) == pytest.approx(19678.1538461538, abs=1e-5)
+    assert float(md["LINE_OFF"]) == pytest.approx(
+        (float(src_md["LINE_OFF"]) - srcyoff + 0.5) * heightratio - 0.5, abs=1e-5
+    )
+    assert float(md["LINE_SCALE"]) == pytest.approx(
+        float(src_md["LINE_SCALE"]) * heightratio, abs=1e-5
+    )
+    assert float(md["SAMP_OFF"]) == pytest.approx(
+        (float(src_md["SAMP_OFF"]) - srcxoff + 0.5) * widthratio - 0.5, abs=1e-5
+    )
+    assert float(md["SAMP_SCALE"]) == pytest.approx(
+        float(src_md["SAMP_SCALE"]) * widthratio, abs=1e-5
     )
 
 

autotest/utilities/test_gdal_translate_lib.py

@@ -1281,11 +1281,15 @@ def test_gdal_translate_ovr_rpc():
     src_rpc = src_ds.GetMetadata("RPC")
     ovr_rpc = ds.GetMetadata("RPC")
     assert ovr_rpc
-    assert float(ovr_rpc["LINE_OFF"]) == pytest.approx(0.5 * float(src_rpc["LINE_OFF"]))
+    assert float(ovr_rpc["LINE_OFF"]) == pytest.approx(
+        0.5 * (float(src_rpc["LINE_OFF"]) + 0.5) - 0.5
+    )
     assert float(ovr_rpc["LINE_SCALE"]) == pytest.approx(
         0.5 * float(src_rpc["LINE_SCALE"])
     )
-    assert float(ovr_rpc["SAMP_OFF"]) == pytest.approx(0.5 * float(src_rpc["SAMP_OFF"]))
+    assert float(ovr_rpc["SAMP_OFF"]) == pytest.approx(
+        0.5 * (float(src_rpc["SAMP_OFF"]) + 0.5) - 0.5
+    )
     assert float(ovr_rpc["SAMP_SCALE"]) == pytest.approx(
         0.5 * float(src_rpc["SAMP_SCALE"])
     )

gcore/gdaloverviewdataset.cpp

@@ -71,7 +71,8 @@ class GDALOverviewDataset final : public GDALDataset
     GDALOverviewBand *m_poMaskBand = nullptr;
 
     static void Rescale(char **&papszMD, const char *pszItem, double dfRatio,
-                        double dfDefaultVal);
+                        double dfDefaultVal, double dfPreShift = 0,
+                        double dfPostShift = 0);
 
   protected:
     CPLErr IRasterIO(GDALRWFlag, int, int, int, int, void *, int, int,
@@ -455,12 +456,17 @@ const GDAL_GCP *GDALOverviewDataset::GetGCPs()
 /*                             Rescale()                                */
 /************************************************************************/
 
+/* static */
 void GDALOverviewDataset::Rescale(char **&papszMD, const char *pszItem,
-                                  double dfRatio, double dfDefaultVal)
+                                  double dfRatio, double dfDefaultVal,
+                                  double dfPreShift /*= 0*/,
+                                  double dfPostShift /*= 0*/)
 {
     double dfVal = CPLAtofM(CSLFetchNameValueDef(
         papszMD, pszItem, CPLSPrintf("%.18g", dfDefaultVal)));
+    dfVal += dfPreShift;
     dfVal *= dfRatio;
+    dfVal += dfPostShift;
     papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.18g", dfVal));
 }
 
@@ -487,18 +493,20 @@ char **GDALOverviewDataset::GetMetadata(const char *pszDomain)
             return papszMD_RPC;
         papszMD_RPC = CSLDuplicate(papszMD);
 
-        Rescale(papszMD_RPC, RPC_LINE_OFF,
-                static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
-                0.0);
-        Rescale(papszMD_RPC, RPC_LINE_SCALE,
-                static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
-                1.0);
-        Rescale(papszMD_RPC, RPC_SAMP_OFF,
-                static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
-                0.0);
-        Rescale(papszMD_RPC, RPC_SAMP_SCALE,
-                static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
-                1.0);
+        const double dfXRatio =
+            static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize();
+        const double dfYRatio =
+            static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize();
+
+        // For line offset and pixel offset, we need to convert from RPC
+        // pixel center registration convention to GDAL pixel top-left corner
+        // registration convention by adding an initial 0.5 shift, and un-apply
+        // it after scaling.
+
+        Rescale(papszMD_RPC, RPC_LINE_OFF, dfYRatio, 0.0, 0.5, -0.5);
+        Rescale(papszMD_RPC, RPC_LINE_SCALE, dfYRatio, 1.0);
+        Rescale(papszMD_RPC, RPC_SAMP_OFF, dfXRatio, 0.0, 0.5, -0.5);
+        Rescale(papszMD_RPC, RPC_SAMP_SCALE, dfXRatio, 1.0);
 
         papszMD = papszMD_RPC;
     }