Merge pull request #10482 from rouault/gti_consistency_index_geometry…

…_vs_source_extent GTI: emit warning if bounding box from the tile's feature extent in the index does not match actual bounding box of tile
OSGeo · Aug 10, 2024 · 4b04099 · 4b04099
2 parents 13ee9c9 + 99ed22e
commit 4b04099
Show file tree

Hide file tree

Showing 2 changed files with 103 additions and 13 deletions.
diff --git a/autotest/gdrivers/gti.py b/autotest/gdrivers/gti.py
@@ -2080,6 +2080,47 @@ def test_gti_overlapping_sources_mask_band(tmp_vsimem):
     ) == (255, 254)
 
 
+def test_gti_consistency_index_geometry_vs_source_extent(tmp_vsimem):
+
+    filename1 = str(tmp_vsimem / "test.tif")
+    ds = gdal.GetDriverByName("GTiff").Create(filename1, 10, 10)
+    ds.SetGeoTransform([2, 1, 0, 49, 0, -1])
+    ds.GetRasterBand(1).Fill(255)
+    expected_cs = ds.GetRasterBand(1).Checksum()
+    del ds
+
+    index_filename = str(tmp_vsimem / "index.gti.gpkg")
+    index_ds, _ = create_basic_tileindex(
+        index_filename,
+        [gdal.Open(filename1)],
+    )
+    del index_ds
+
+    vrt_ds = gdal.Open(index_filename)
+    with gdal.quiet_errors():
+        gdal.ErrorReset()
+        assert vrt_ds.GetRasterBand(1).Checksum() == expected_cs
+        assert gdal.GetLastErrorMsg() == ""
+
+    # No intersection
+    with gdal.Open(filename1, gdal.GA_Update) as ds:
+        ds.SetGeoTransform([100, 1, 0, 49, 0, -1])
+
+    vrt_ds = gdal.Open(index_filename)
+    with gdal.quiet_errors():
+        assert vrt_ds.GetRasterBand(1).Checksum() == 0
+        assert "does not intersect at all" in gdal.GetLastErrorMsg()
+
+    # Partial intersection
+    with gdal.Open(filename1, gdal.GA_Update) as ds:
+        ds.SetGeoTransform([4, 1, 0, 49, 0, -1])
+
+    vrt_ds = gdal.Open(index_filename)
+    with gdal.quiet_errors():
+        assert vrt_ds.GetRasterBand(1).Checksum() == 958
+        assert "does not fully contain" in gdal.GetLastErrorMsg()
+
+
 def test_gti_mask_band_explicit(tmp_vsimem):
 
     index_filename = str(tmp_vsimem / "index.gti.gpkg")

diff --git a/frmts/vrt/gdaltileindexdataset.cpp b/frmts/vrt/gdaltileindexdataset.cpp
@@ -3174,19 +3174,13 @@ bool GDALTileIndexDataset::GetSourceDesc(const std::string &osTileName,
         poSource = std::move(poComplexSource);
     }
 
-    if (!GetSrcDstWin(adfGeoTransformTile, poTileDS->GetRasterXSize(),
-                      poTileDS->GetRasterYSize(), m_adfGeoTransform.data(),
-                      GetRasterXSize(), GetRasterYSize(),
-                      &poSource->m_dfSrcXOff, &poSource->m_dfSrcYOff,
-                      &poSource->m_dfSrcXSize, &poSource->m_dfSrcYSize,
-                      &poSource->m_dfDstXOff, &poSource->m_dfDstYOff,
-                      &poSource->m_dfDstXSize, &poSource->m_dfDstYSize))
-    {
-        // Should not happen on a consistent tile index
-        CPLDebug("VRT", "Tile %s does not actually intersect area of interest",
-                 osTileName.c_str());
-        return false;
-    }
+    GetSrcDstWin(adfGeoTransformTile, poTileDS->GetRasterXSize(),
+                 poTileDS->GetRasterYSize(), m_adfGeoTransform.data(),
+                 GetRasterXSize(), GetRasterYSize(), &poSource->m_dfSrcXOff,
+                 &poSource->m_dfSrcYOff, &poSource->m_dfSrcXSize,
+                 &poSource->m_dfSrcYSize, &poSource->m_dfDstXOff,
+                 &poSource->m_dfDstYOff, &poSource->m_dfDstXSize,
+                 &poSource->m_dfDstYSize);
 
     oSourceDesc.osName = osTileName;
     oSourceDesc.poDS = std::move(poTileDS);
@@ -3390,6 +3384,61 @@ bool GDALTileIndexDataset::CollectSources(double dfXOff, double dfYOff,
         if (!GetSourceDesc(osTileName, oSourceDesc, nullptr))
             continue;
 
+        // Check consistency of bounding box in tile index vs actual
+        // extent of the tile.
+        double adfTileGT[6];
+        if (oSourceDesc.poDS->GetGeoTransform(adfTileGT) == CE_None &&
+            adfTileGT[GT_ROTATION_PARAM1] == 0 &&
+            adfTileGT[GT_ROTATION_PARAM2] == 0)
+        {
+            OGREnvelope sActualTileExtent;
+            sActualTileExtent.MinX = adfTileGT[GT_TOPLEFT_X];
+            sActualTileExtent.MaxX =
+                sActualTileExtent.MinX +
+                oSourceDesc.poDS->GetRasterXSize() * adfTileGT[GT_WE_RES];
+            sActualTileExtent.MaxY = adfTileGT[GT_TOPLEFT_Y];
+            sActualTileExtent.MinY =
+                sActualTileExtent.MaxY +
+                oSourceDesc.poDS->GetRasterYSize() * adfTileGT[GT_NS_RES];
+            const auto poGeom = poFeature->GetGeometryRef();
+            if (poGeom && !poGeom->IsEmpty())
+            {
+                OGREnvelope sGeomTileExtent;
+                poGeom->getEnvelope(&sGeomTileExtent);
+                sGeomTileExtent.MinX -= m_adfGeoTransform[GT_WE_RES];
+                sGeomTileExtent.MaxX += m_adfGeoTransform[GT_WE_RES];
+                sGeomTileExtent.MinY -= std::fabs(m_adfGeoTransform[GT_NS_RES]);
+                sGeomTileExtent.MaxY += std::fabs(m_adfGeoTransform[GT_NS_RES]);
+                if (!sGeomTileExtent.Contains(sActualTileExtent))
+                {
+                    if (!sGeomTileExtent.Intersects(sActualTileExtent))
+                    {
+                        CPLError(CE_Warning, CPLE_AppDefined,
+                                 "Tile index is out of sync with actual "
+                                 "extent of %s. Bounding box from tile index "
+                                 "is (%g, %g, %g, %g) does not intersect at "
+                                 "all bounding box from tile (%g, %g, %g, %g)",
+                                 osTileName.c_str(), sGeomTileExtent.MinX,
+                                 sGeomTileExtent.MinY, sGeomTileExtent.MaxX,
+                                 sGeomTileExtent.MaxY, sActualTileExtent.MinX,
+                                 sActualTileExtent.MinY, sActualTileExtent.MaxX,
+                                 sActualTileExtent.MaxY);
+                        continue;
+                    }
+                    CPLError(CE_Warning, CPLE_AppDefined,
+                             "Tile index is out of sync with actual extent "
+                             "of %s. Bounding box from tile index is (%g, %g, "
+                             "%g, %g) does not fully contain bounding box from "
+                             "tile (%g, %g, %g, %g)",
+                             osTileName.c_str(), sGeomTileExtent.MinX,
+                             sGeomTileExtent.MinY, sGeomTileExtent.MaxX,
+                             sGeomTileExtent.MaxY, sActualTileExtent.MinX,
+                             sActualTileExtent.MinY, sActualTileExtent.MaxX,
+                             sActualTileExtent.MaxY);
+                }
+            }
+        }
+
         const auto &poSource = oSourceDesc.poSource;
         if (dfXOff >= poSource->m_dfDstXOff + poSource->m_dfDstXSize ||
             dfYOff >= poSource->m_dfDstYOff + poSource->m_dfDstYSize ||