Merge pull request OSGeo#12814 from dbaston/raster-clip-fix-bbox-check

gdal raster clip: fix bbox check

Author: rouault

apps/gdalalg_raster_clip.cpp

@@ -252,11 +252,9 @@ bool GDALRasterClipAlgorithm::RunStep(GDALPipelineStepRunContext &)
             poClipGeomInSrcSRS->getEnvelope(&env);
         }
 
-        if (!m_allowExtentOutsideSource &&
-            !(env.MinX >= gt[0] &&
-              env.MaxX <= gt[0] + gt[1] * poSrcDS->GetRasterXSize() &&
-              env.MaxY >= gt[3] &&
-              env.MinY <= gt[3] + gt[5] * poSrcDS->GetRasterYSize()))
+        OGREnvelope rasterEnv;
+        poSrcDS->GetExtent(&rasterEnv, nullptr);
+        if (!m_allowExtentOutsideSource && !rasterEnv.Contains(env))
         {
             ReportError(CE_Failure, CPLE_AppDefined,
                         "Clipping geometry is partially or totally outside the "

autotest/gcore/basic_test.py

@@ -19,7 +19,7 @@
 import gdaltest
 import pytest
 
-from osgeo import gdal
+from osgeo import gdal, osr
 
 
 ###############################################################################
@@ -1108,3 +1108,56 @@ def test_basic_block_windows(tmp_vsimem):
         assert windows[8].yoff == 512
         assert windows[8].xsize == 1050 - 1024
         assert windows[8].ysize == 600 - 512
+
+
+###############################################################################
+# Test GetExtent()
+
+
+def test_basic_get_extent():
+
+    with gdal.Open("data/byte.tif") as ds:
+        assert ds.GetExtent() == (440720.0, 441920.0, 3750120.0, 3751320.0)
+
+
+def test_basic_get_extent_reprojected():
+
+    wgs84 = osr.SpatialReference(epsg=4326)
+    wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+
+    with gdal.Open("data/byte.tif") as ds:
+        assert ds.GetExtent(wgs84) == pytest.approx(
+            (-117.642, -117.629, 33.892, 33.902), abs=1e-3
+        )
+
+
+def test_basic_get_extent_no_crs():
+
+    with gdal.GetDriverByName("MEM").Create("", 5, 5) as ds:
+        ds.SetGeoTransform((3, 0.5, 0, 7, 0, -1))
+        assert ds.GetExtent() == (3, 5.5, 2, 7)
+
+
+def test_basic_get_extent_no_crs_reprojected():
+
+    wgs84 = osr.SpatialReference(epsg=4326)
+    wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+
+    with gdal.GetDriverByName("MEM").Create("", 5, 5) as ds:
+        ds.SetGeoTransform((3, 0.5, 0, 7, 0, -1))
+        assert ds.GetExtent(wgs84) is None
+
+
+def test_basic_get_extent_bottom_up():
+
+    with gdal.GetDriverByName("MEM").Create("", 5, 5) as ds:
+        ds.SetGeoTransform((3, 0.5, 0, 7, 0, 1))
+        assert ds.GetExtent() == (3, 5.5, 7, 12)
+
+
+def test_basic_get_extent_rotated():
+
+    with gdal.Open("data/geomatrix.tif") as ds:
+        assert ds.GetExtent() == pytest.approx(
+            (1840900, 1841030, 1143870, 1144000), abs=4
+        )

autotest/utilities/test_gdalalg_raster_clip.py

@@ -169,6 +169,28 @@ def test_gdalalg_raster_clip_bbox_crs():
     )
 
 
+def test_gdalgalg_raster_clip_geometry(tmp_vsimem):
+
+    alg = get_alg()
+
+    alg["input"] = "../gcore/data/byte.tif"
+    alg["output"] = ""
+    alg["output-format"] = "MEM"
+    alg[
+        "geometry"
+    ] = "POLYGON ((440885 3750741, 441344 3750294, 441612 3750501, 441773 3751203, 441545 3751254, 441576 3750847, 441576 3750847, 440885 3750741))"
+
+    assert alg.Run()
+
+    ds = alg["output"].GetDataset()
+    assert ds.RasterXSize == 16
+    assert ds.RasterYSize == 17
+    assert ds.GetSpatialRef().GetAuthorityCode(None) == "26711"
+    assert ds.GetGeoTransform() == pytest.approx(
+        (440840, 60.0, 0.0, 3751260, 0.0, -60.0), rel=1e-8
+    )
+
+
 def test_gdalalg_raster_clip_geometry_add_alpha():
 
     src_ds = gdal.Open("../gcore/data/byte.tif")

gcore/gdaldataset.cpp

@@ -10621,7 +10621,7 @@ CPLErr GDALDataset::GetExtent(OGREnvelope *psExtent,
             if (auto poLayer = poThisDS->GetLayer(0))
                 poThisCRS = poLayer->GetSpatialRef();
         }
-        if (!poThisCRS)
+        if (poCRS && !poThisCRS)
             return CE_Failure;
     }
     if (!poCRS)
@@ -10633,34 +10633,46 @@ CPLErr GDALDataset::GetExtent(OGREnvelope *psExtent,
     const bool bHasGT = poThisDS->GetGeoTransform(gt) == CE_None;
     if (bHasGT)
     {
-        auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
-            OGRCreateCoordinateTransformation(poThisCRS, poCRS));
+        std::unique_ptr<OGRCoordinateTransformation> poCT;
+        if (poCRS)
+        {
+            poCT.reset(OGRCreateCoordinateTransformation(poThisCRS, poCRS));
+        }
+
+        constexpr int DENSIFY_POINT_COUNT = 21;
+        double dfULX = gt[0];
+        double dfULY = gt[3];
+        double dfURX = 0, dfURY = 0;
+        gt.Apply(nRasterXSize, 0, &dfURX, &dfURY);
+        double dfLLX = 0, dfLLY = 0;
+        gt.Apply(0, nRasterYSize, &dfLLX, &dfLLY);
+        double dfLRX = 0, dfLRY = 0;
+        gt.Apply(nRasterXSize, nRasterYSize, &dfLRX, &dfLRY);
+        const double xmin =
+            std::min(std::min(dfULX, dfURX), std::min(dfLLX, dfLRX));
+        const double ymin =
+            std::min(std::min(dfULY, dfURY), std::min(dfLLY, dfLRY));
+        const double xmax =
+            std::max(std::max(dfULX, dfURX), std::max(dfLLX, dfLRX));
+        const double ymax =
+            std::max(std::max(dfULY, dfURY), std::max(dfLLY, dfLRY));
         if (poCT)
         {
-            constexpr int DENSIFY_POINT_COUNT = 21;
             OGREnvelope sEnvTmp;
-            double dfULX = gt[0];
-            double dfULY = gt[3];
-            double dfURX = 0, dfURY = 0;
-            gt.Apply(nRasterXSize, 0, &dfURX, &dfURY);
-            double dfLLX = 0, dfLLY = 0;
-            gt.Apply(0, nRasterYSize, &dfLLX, &dfLLY);
-            double dfLRX = 0, dfLRY = 0;
-            gt.Apply(nRasterXSize, nRasterYSize, &dfLRX, &dfLRY);
-            const double xmin =
-                std::min(std::min(dfULX, dfURX), std::min(dfLLX, dfLRX));
-            const double ymin =
-                std::min(std::min(dfULY, dfURY), std::min(dfLLY, dfLRY));
-            const double xmax =
-                std::max(std::max(dfULX, dfURX), std::max(dfLLX, dfLRX));
-            const double ymax =
-                std::max(std::max(dfULY, dfURY), std::max(dfLLY, dfLRY));
-            if (poCT->TransformBounds(xmin, ymin, xmax, ymax, &(sEnvTmp.MinX),
-                                      &(sEnvTmp.MinY), &(sEnvTmp.MaxX),
-                                      &(sEnvTmp.MaxY), DENSIFY_POINT_COUNT))
+            if (!poCT->TransformBounds(xmin, ymin, xmax, ymax, &(sEnvTmp.MinX),
+                                       &(sEnvTmp.MinY), &(sEnvTmp.MaxX),
+                                       &(sEnvTmp.MaxY), DENSIFY_POINT_COUNT))
             {
-                *psExtent = std::move(sEnvTmp);
+                return CE_Failure;
             }
+            *psExtent = sEnvTmp;
+        }
+        else
+        {
+            psExtent->MinX = xmin;
+            psExtent->MinY = ymin;
+            psExtent->MaxX = xmax;
+            psExtent->MaxY = ymax;
         }
     }