Merge pull request #13418 from rouault/fix_13416

GTI: accept tiles with south-up orientation (auto-warp them to north-up)

Author: rouault

autotest/gdrivers/gti.py

@@ -329,7 +329,7 @@ def test_gti_prototype_tile_wrong_gt_5th_value(tmp_vsimem):
         gdal.Open(index_filename)
 
 
-def test_gti_prototype_tile_wrong_gt_6th_value(tmp_vsimem):
+def test_gti_prototype_tile_wrong_gt_2nd_value(tmp_vsimem):
 
     index_filename = str(tmp_vsimem / "index.gti.gpkg")
     protods_filename = str(tmp_vsimem / "protods_filename.tif")
@@ -346,7 +346,7 @@ def test_gti_prototype_tile_wrong_gt_6th_value(tmp_vsimem):
     lyr.CreateFeature(f)
     del index_ds
 
-    with pytest.raises(Exception, match="6th value of GeoTransform"):
+    with pytest.raises(Exception, match="2nd value of GeoTransform"):
         gdal.Open(index_filename)
 
 
@@ -369,6 +369,28 @@ def test_gti_no_extent(tmp_vsimem):
         gdal.Open(index_filename)
 
 
+def test_gti_south_up(tmp_vsimem):
+
+    index_filename = str(tmp_vsimem / "index.gti.gpkg")
+    protods_filename = str(tmp_vsimem / "protods_filename.tif")
+    with gdal.GetDriverByName("GTiff").Create(protods_filename, 1, 2) as ds:
+        ds.SetGeoTransform([0, 1, 0, 10, 0, 1])
+        ds.WriteRaster(0, 0, 1, 2, b"\x01\x02")
+
+    index_ds = ogr.GetDriverByName("GPKG").CreateDataSource(index_filename)
+    lyr = index_ds.CreateLayer("index")
+    lyr.CreateField(ogr.FieldDefn("location"))
+    f = ogr.Feature(lyr.GetLayerDefn())
+    f["location"] = protods_filename
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON((0 10,0 12,1 12,1 12,0 12))"))
+    lyr.CreateFeature(f)
+    del index_ds
+
+    ds = gdal.Open(index_filename)
+    assert ds.GetGeoTransform() == (0, 1, 0, 12, 0, -1)
+    assert ds.ReadRaster() == b"\x02\x01"
+
+
 def test_gti_too_big_x(tmp_vsimem):
 
     index_filename = str(tmp_vsimem / "index.gti.gpkg")

frmts/gti/gdaltileindexdataset.cpp

@@ -1553,16 +1553,23 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
                      "5th value of GeoTransform of %s must be 0", pszTileName);
             return false;
         }
-        if (!(gtTile[GT_NS_RES] < 0))
+
+        const double dfResX = gtTile[GT_WE_RES];
+        const double dfResY = gtTile[GT_NS_RES];
+        if (!(dfResX > 0))
         {
             CPLError(CE_Failure, CPLE_AppDefined,
-                     "6th value of GeoTransform of %s must be < 0",
+                     "2nd value of GeoTransform of %s must be > 0",
+                     pszTileName);
+            return false;
+        }
+        if (!(dfResY != 0))
+        {
+            CPLError(CE_Failure, CPLE_AppDefined,
+                     "6th value of GeoTransform of %s must be != 0",
                      pszTileName);
             return false;
         }
-
-        const double dfResX = gtTile[GT_WE_RES];
-        const double dfResY = -gtTile[GT_NS_RES];
 
         if (!sEnvelope.IsInit() &&
             m_poLayer->GetExtent(&sEnvelope, /* bForce = */ false) ==
@@ -1587,7 +1594,8 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
             return false;
         }
 
-        const double dfYSize = (sEnvelope.MaxY - sEnvelope.MinY) / dfResY;
+        const double dfYSize =
+            (sEnvelope.MaxY - sEnvelope.MinY) / std::fabs(dfResY);
         if (!(dfYSize >= 0 && dfYSize < INT_MAX))
         {
             CPLError(CE_Failure, CPLE_AppDefined,
@@ -1600,7 +1608,8 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
         m_gt[GT_ROTATION_PARAM1] = 0;
         m_gt[GT_TOPLEFT_Y] = sEnvelope.MaxY;
         m_gt[GT_ROTATION_PARAM2] = 0;
-        m_gt[GT_NS_RES] = -dfResY;
+        m_gt[GT_NS_RES] = -std::fabs(dfResY);
+
         nRasterXSize = static_cast<int>(std::ceil(dfXSize));
         nRasterYSize = static_cast<int>(std::ceil(dfYSize));
     }
@@ -1637,6 +1646,12 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
         {
             m_gt[i] = CPLAtof(aosTokens[i]);
         }
+        if (!(m_gt[GT_WE_RES] > 0))
+        {
+            CPLError(CE_Failure, CPLE_AppDefined, "2nd value of %s must be > 0",
+                     MD_GEOTRANSFORM);
+            return false;
+        }
         if (!(m_gt[GT_ROTATION_PARAM1] == 0))
         {
             CPLError(CE_Failure, CPLE_AppDefined, "3rd value of %s must be 0",
@@ -1655,7 +1670,6 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
                      MD_GEOTRANSFORM);
             return false;
         }
-
         nRasterXSize = nXSize;
         nRasterYSize = nYSize;
     }
@@ -3492,16 +3506,37 @@ bool GDALTileIndexDataset::GetSourceDesc(const std::string &osTileName,
         if (!GTIDoPaletteExpansionIfNeeded(poTileDS, nBands))
             return false;
 
-        const OGRSpatialReference *poTileSRS;
-        if (!m_oSRS.IsEmpty() &&
-            (poTileSRS = poTileDS->GetSpatialRef()) != nullptr &&
+        bool bWarpVRT = false;
+        bool bExportSRS = false;
+        bool bAddAlphaToVRT = false;
+        const OGRSpatialReference *poTileSRS = poTileDS->GetSpatialRef();
+        GDALGeoTransform tileGT;
+        if (!m_oSRS.IsEmpty() && poTileSRS != nullptr &&
             !m_oSRS.IsSame(poTileSRS))
         {
             CPLDebug("VRT",
                      "Tile %s has not the same SRS as the VRT. "
                      "Proceed to on-the-fly warping",
                      osTileName.c_str());
+            bWarpVRT = true;
+            bExportSRS = true;
+            bAddAlphaToVRT = true;
+        }
+        else if (poTileDS->GetGeoTransform(tileGT) == CE_None &&
+                 tileGT[GT_NS_RES] > 0 &&
+                 ((m_oSRS.IsEmpty() && poTileSRS == nullptr) ||
+                  (!m_oSRS.IsEmpty() && poTileSRS && m_oSRS.IsSame(poTileSRS))))
 
+        {
+            CPLDebug("VRT",
+                     "Tile %s is south-up oriented. "
+                     "Proceed to on-the-fly warping",
+                     osTileName.c_str());
+            bWarpVRT = true;
+        }
+
+        if (bWarpVRT)
+        {
             CPLStringList aosOptions;
             aosOptions.AddString("-of");
             aosOptions.AddString("VRT");
@@ -3514,25 +3549,29 @@ bool GDALTileIndexDataset::GetSourceDesc(const std::string &osTileName,
                 aosOptions.AddString(m_osResampling.c_str());
             }
 
-            if (m_osWKT.empty())
+            if (bExportSRS)
             {
-                char *pszWKT = nullptr;
-                const char *const apszWKTOptions[] = {"FORMAT=WKT2_2019",
-                                                      nullptr};
-                m_oSRS.exportToWkt(&pszWKT, apszWKTOptions);
-                if (pszWKT)
-                    m_osWKT = pszWKT;
-                CPLFree(pszWKT);
-            }
-            if (m_osWKT.empty())
-            {
-                CPLError(CE_Failure, CPLE_AppDefined,
-                         "Cannot export VRT SRS to WKT2");
-                return false;
-            }
+                if (m_osWKT.empty())
+                {
+                    char *pszWKT = nullptr;
+                    const char *const apszWKTOptions[] = {"FORMAT=WKT2_2019",
+                                                          nullptr};
+                    m_oSRS.exportToWkt(&pszWKT, apszWKTOptions);
+                    if (pszWKT)
+                        m_osWKT = pszWKT;
+                    CPLFree(pszWKT);
+
+                    if (m_osWKT.empty())
+                    {
+                        CPLError(CE_Failure, CPLE_AppDefined,
+                                 "Cannot export VRT SRS to WKT2");
+                        return false;
+                    }
+                }
 
-            aosOptions.AddString("-t_srs");
-            aosOptions.AddString(m_osWKT.c_str());
+                aosOptions.AddString("-t_srs");
+                aosOptions.AddString(m_osWKT.c_str());
+            }
 
             // First pass to get the extent of the tile in the
             // target VRT SRS
@@ -3592,7 +3631,8 @@ bool GDALTileIndexDataset::GetSourceDesc(const std::string &osTileName,
             aosOptions.AddString(CPLSPrintf("%.17g", dfVRTResX));
             aosOptions.AddString(CPLSPrintf("%.17g", dfVRTResYAbs));
 
-            aosOptions.AddString("-dstalpha");
+            if (bAddAlphaToVRT)
+                aosOptions.AddString("-dstalpha");
 
             psWarpOptions = GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
             poWarpDS.reset(GDALDataset::FromHandle(GDALWarp(

gcore/gdalproxypool.cpp

@@ -1110,8 +1110,11 @@ const OGRSpatialReference *GDALProxyPoolDataset::GetGCPSpatialRef() const
     if (poUnderlyingDataset == nullptr)
         return nullptr;
 
-    m_poGCPSRS->Release();
-    m_poGCPSRS = nullptr;
+    if (m_poGCPSRS)
+    {
+        m_poGCPSRS->Release();
+        m_poGCPSRS = nullptr;
+    }
 
     const auto poUnderlyingGCPSRS = poUnderlyingDataset->GetGCPSpatialRef();
     if (poUnderlyingGCPSRS)