GDALDataset::GetExtent: Return extent for undefined CRS

Author: dbaston

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
+        )

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;
         }
     }