Merge pull request #8436 from rouault/gdalwarp_ct_cutline

gdalwarp: fix error when using -ct and -cutline (master only), and actually use the -ct when possible if sourceCRS != targetCRS and targetCRS == cutlineCRS
OSGeo · Sep 20, 2023 · 8122e2c · 8122e2c
2 parents d8086cb + 93f77d7
commit 8122e2c
Show file tree

Hide file tree

Showing 2 changed files with 107 additions and 18 deletions.
diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp
@@ -400,6 +400,40 @@ static CPLString GetSrcDSProjection(GDALDatasetH hDS, CSLConstList papszTO)
     return pszProjection ? pszProjection : "";
 }
 
+/************************************************************************/
+/*                      CreateCTCutlineToSrc()                          */
+/************************************************************************/
+
+static std::unique_ptr<OGRCoordinateTransformation> CreateCTCutlineToSrc(
+    const OGRSpatialReference *poRasterSRS, const OGRSpatialReference *poDstSRS,
+    const OGRSpatialReference *poCutlineSRS, CSLConstList papszTO)
+{
+    const OGRSpatialReference *poCutlineOrTargetSRS =
+        poCutlineSRS ? poCutlineSRS : poDstSRS;
+    std::unique_ptr<OGRCoordinateTransformation> poCTCutlineToSrc;
+    if (poCutlineOrTargetSRS && poRasterSRS &&
+        !poCutlineOrTargetSRS->IsSame(poRasterSRS))
+    {
+        OGRCoordinateTransformationOptions oOptions;
+        // If the cutline SRS is the same as the target SRS and there is
+        // an explicit -ct between the source SRS and the target SRS, then
+        // use it in the reverse way to transform from the cutline SRS to
+        // the source SRS.
+        if (poDstSRS && poCutlineOrTargetSRS->IsSame(poDstSRS))
+        {
+            const char *pszCT =
+                CSLFetchNameValue(papszTO, "COORDINATE_OPERATION");
+            if (pszCT)
+            {
+                oOptions.SetCoordinateOperation(pszCT, /* bInverse = */ true);
+            }
+        }
+        poCTCutlineToSrc.reset(OGRCreateCoordinateTransformation(
+            poCutlineOrTargetSRS, poRasterSRS, oOptions));
+    }
+    return poCTCutlineToSrc;
+}
+
 /************************************************************************/
 /*                           CropToCutline()                            */
 /************************************************************************/
@@ -470,16 +504,10 @@ static CPLErr CropToCutline(const OGRGeometry *poCutline, CSLConstList papszTO,
     }
 
     auto poCutlineGeom = std::unique_ptr<OGRGeometry>(poCutline->clone());
-    const OGRSpatialReference *poCutlineOrTargetSRS =
-        poCutlineSRS ? poCutlineSRS : poDstSRS.get();
-    std::unique_ptr<OGRCoordinateTransformation> poCTCutlineToSrc;
-    std::unique_ptr<OGRCoordinateTransformation> poCTSrcToDst;
+    auto poCTCutlineToSrc = CreateCTCutlineToSrc(poSrcSRS.get(), poDstSRS.get(),
+                                                 poCutlineSRS, papszTO);
 
-    if (!poCutlineOrTargetSRS->IsSame(poSrcSRS.get()))
-    {
-        poCTCutlineToSrc.reset(OGRCreateCoordinateTransformation(
-            poCutlineOrTargetSRS, poSrcSRS.get()));
-    }
+    std::unique_ptr<OGRCoordinateTransformation> poCTSrcToDst;
     if (!poSrcSRS->IsSame(poDstSRS.get()))
     {
         poCTSrcToDst.reset(
@@ -4775,15 +4803,8 @@ static CPLErr TransformCutlineToSource(GDALDataset *poSrcDS,
                  "Cutline results may be incorrect.");
     }
 
-    const OGRSpatialReference *poCutlineOrTargetSRS =
-        poCutlineSRS ? poCutlineSRS : poDstSRS.get();
-    std::unique_ptr<OGRCoordinateTransformation> poCTCutlineToSrc;
-    if (poCutlineOrTargetSRS && poRasterSRS &&
-        !poCutlineOrTargetSRS->IsSame(poRasterSRS.get()))
-    {
-        poCTCutlineToSrc.reset(OGRCreateCoordinateTransformation(
-            poCutlineOrTargetSRS, poRasterSRS.get()));
-    }
+    auto poCTCutlineToSrc = CreateCTCutlineToSrc(
+        poRasterSRS.get(), poDstSRS.get(), poCutlineSRS, papszTO_In);
 
     CPLStringList aosTO(papszTO_In);
 
@@ -4792,6 +4813,9 @@ static CPLErr TransformCutlineToSource(GDALDataset *poSrcDS,
         // Avoid any reprojection when using the GenImgProjTransformer
         aosTO.SetNameValue("DST_SRS", osProjection.c_str());
     }
+    aosTO.SetNameValue("SRC_COORDINATE_EPOCH", nullptr);
+    aosTO.SetNameValue("DST_COORDINATE_EPOCH", nullptr);
+    aosTO.SetNameValue("COORDINATE_OPERATION", nullptr);
 
     /* -------------------------------------------------------------------- */
     /*      It may be unwise to let the mask geometry be re-wrapped by      */
@@ -4834,7 +4858,12 @@ static CPLErr TransformCutlineToSource(GDALDataset *poSrcDS,
         }
     }
     if (poMultiPolygon->transform(&oTransformer) != OGRERR_NONE)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "poMultiPolygon->transform(&oTransformer) failed at line %d",
+                 __LINE__);
         eErr = OGRERR_FAILURE;
+    }
     const double dfInitialMaxLengthInPixels =
         GetMaximumSegmentLength(poMultiPolygon.get());
 

diff --git a/autotest/utilities/test_gdalwarp_lib.py b/autotest/utilities/test_gdalwarp_lib.py
@@ -472,6 +472,66 @@ def test_gdalwarp_lib_21():
     ds = None
 
 
+###############################################################################
+# Test cutline with sourceCRS != targetCRS and targetCRS == cutlineCRS
+
+
+@pytest.mark.require_driver("CSV")
+@pytest.mark.require_driver("GPKG")
+def test_gdalwarp_lib_cutline_reprojection(tmp_vsimem):
+
+    cutline_filename = str(tmp_vsimem / "cutline.gpkg")
+    gdal.VectorTranslate(
+        cutline_filename, "data/cutline.vrt", dstSRS="EPSG:4267", reproject=True
+    )
+
+    ds = gdal.Warp(
+        "",
+        "../gcore/data/utmsmall.tif",
+        format="MEM",
+        dstSRS="EPSG:4267",
+        cutlineDSName=cutline_filename,
+    )
+    assert ds is not None
+
+    assert ds.GetRasterBand(1).Checksum() == 18883, "Bad checksum"
+
+    ds = None
+
+
+###############################################################################
+# Test cutline with sourceCRS != targetCRS and targetCRS == cutlineCRS and coordinateOperation
+
+
+@pytest.mark.require_driver("CSV")
+@pytest.mark.require_driver("GPKG")
+def test_gdalwarp_lib_cutline_reprojection_and_coordinate_operation(tmp_vsimem):
+
+    cutline_filename = str(tmp_vsimem / "cutline.gpkg")
+    ct = "+proj=pipeline +step +proj=affine +xoff=10000 +step +inv +proj=utm +zone=11 +ellps=clrk66 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"
+    gdal.VectorTranslate(
+        cutline_filename,
+        "data/cutline.vrt",
+        dstSRS="EPSG:4267",
+        coordinateOperation=ct,
+        reproject=True,
+    )
+
+    ds = gdal.Warp(
+        "",
+        "../gcore/data/utmsmall.tif",
+        format="MEM",
+        dstSRS="EPSG:4267",
+        cutlineDSName=cutline_filename,
+        coordinateOperation=ct,
+    )
+    assert ds is not None
+
+    assert ds.GetRasterBand(1).Checksum() == 18914, "Bad checksum"
+
+    ds = None
+
+
 ###############################################################################
 # Test cutline from PostGIS (mostly to check that we open the dataset in
 # vector mode, and not in raster mode, which would cause the PostGISRaster