Merge pull request #8847 from OSGeo/backport-8839-to-release/3.8

[Backport release/3.8] gdal_footprint: really fix datasets with alpha band, by vectorizing only the alpha band and not the other ones (really fixes #8792, fixes #8834)
OSGeo · Nov 28, 2023 · cb40027 · cb40027
2 parents 06c629c + 3c51645
commit cb40027
Show file tree

Hide file tree

Showing 3 changed files with 91 additions and 5 deletions.
diff --git a/apps/gdal_footprint_lib.cpp b/apps/gdal_footprint_lib.cpp
@@ -703,14 +703,13 @@ static bool GDALFootprintProcess(GDALDataset *poSrcDS, OGRLayer *poDstLayer,
         {
             GDALRasterBand *poMaskBand;
             const int nMaskFlags = poBand->GetMaskFlags();
-            if ((nMaskFlags & GMF_ALPHA) != 0 ||
-                poBand->GetColorInterpretation() == GCI_AlphaBand)
+            if (poBand->GetColorInterpretation() == GCI_AlphaBand)
             {
                 poMaskBand = poBand;
             }
             else
             {
-                if ((nMaskFlags & GMF_PER_DATASET) != 0)
+                if ((nMaskFlags & GMF_PER_DATASET) == 0)
                 {
                     bGlobalMask = false;
                 }

diff --git a/autotest/utilities/test_gdal_footprint_lib.py b/autotest/utilities/test_gdal_footprint_lib.py
@@ -439,8 +439,7 @@ def test_gdaldem_footprint_rgba_overviews():
     for i in range(4):
         src_ds.GetRasterBand(i + 1).SetColorInterpretation(gdal.GCI_RedBand + i)
     src_ds.BuildOverviews("NONE", [2])
-    for i in range(4):
-        src_ds.GetRasterBand(i + 1).GetOverview(0).WriteRaster(1, 1, 1, 1, b"\xFF")
+    src_ds.GetRasterBand(4).GetOverview(0).WriteRaster(1, 1, 1, 1, b"\xFF")
     out_ds = gdal.Footprint(
         "",
         src_ds,
@@ -452,3 +451,86 @@ def test_gdaldem_footprint_rgba_overviews():
     lyr = out_ds.GetLayer(0)
     f = lyr.GetNextFeature()
     ogrtest.check_feature_geometry(f, "MULTIPOLYGON (((2 2,2 4,4 4,4 2,2 2)))")
+
+
+###############################################################################
+def test_gdal_footprint_lib_union():
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 1, 2)
+    src_ds.GetRasterBand(1).SetNoDataValue(0)
+    src_ds.GetRasterBand(1).WriteRaster(0, 0, 1, 1, b"\xFF")
+    src_ds.GetRasterBand(2).SetNoDataValue(0)
+    src_ds.GetRasterBand(2).WriteRaster(1, 0, 1, 1, b"\xFF")
+    out_ds = gdal.Footprint(
+        "",
+        src_ds,
+        format="Memory",
+        targetCoordinateSystem="pixel",
+    )
+    assert out_ds is not None
+    lyr = out_ds.GetLayer(0)
+    f = lyr.GetNextFeature()
+    ogrtest.check_feature_geometry(f, "MULTIPOLYGON (((0 0,0 1,2 1,2 0,0 0)))")
+
+
+###############################################################################
+def test_gdal_footprint_lib_intersection_none():
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 1, 2)
+    src_ds.GetRasterBand(1).SetNoDataValue(0)
+    src_ds.GetRasterBand(1).WriteRaster(0, 0, 1, 1, b"\xFF")
+    src_ds.GetRasterBand(2).SetNoDataValue(0)
+    src_ds.GetRasterBand(2).WriteRaster(1, 0, 1, 1, b"\xFF")
+    out_ds = gdal.Footprint(
+        "",
+        src_ds,
+        format="Memory",
+        targetCoordinateSystem="pixel",
+        combineBands="intersection",
+    )
+    assert out_ds is not None
+    lyr = out_ds.GetLayer(0)
+    f = lyr.GetNextFeature()
+    assert f is None
+
+
+###############################################################################
+def test_gdal_footprint_lib_intersection_partial():
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 1, 2)
+    src_ds.GetRasterBand(1).SetNoDataValue(0)
+    src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, b"\xFF\xFF")
+    src_ds.GetRasterBand(2).SetNoDataValue(0)
+    src_ds.GetRasterBand(2).WriteRaster(1, 0, 1, 1, b"\xFF")
+    out_ds = gdal.Footprint(
+        "",
+        src_ds,
+        format="Memory",
+        targetCoordinateSystem="pixel",
+        combineBands="intersection",
+    )
+    assert out_ds is not None
+    lyr = out_ds.GetLayer(0)
+    f = lyr.GetNextFeature()
+    ogrtest.check_feature_geometry(f, "MULTIPOLYGON (((1 0,1 1,2 1,2 0,1 0)))")
+
+
+###############################################################################
+def test_gdal_footprint_lib_intersection_full():
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 1, 2)
+    src_ds.GetRasterBand(1).SetNoDataValue(0)
+    src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, b"\xFF\xFF")
+    src_ds.GetRasterBand(2).SetNoDataValue(0)
+    src_ds.GetRasterBand(2).WriteRaster(0, 0, 2, 1, b"\xFF\xFF")
+    out_ds = gdal.Footprint(
+        "",
+        src_ds,
+        format="Memory",
+        targetCoordinateSystem="pixel",
+        combineBands="intersection",
+    )
+    assert out_ds is not None
+    lyr = out_ds.GetLayer(0)
+    f = lyr.GetNextFeature()
+    ogrtest.check_feature_geometry(f, "MULTIPOLYGON (((0 0,0 1,2 1,2 0,0 0)))")
diff --git a/swig/include/python/gdal_python.i b/swig/include/python/gdal_python.i
@@ -3406,6 +3406,7 @@ def Rasterize(destNameOrDestDS, srcDS, **kwargs):
 def FootprintOptions(options=None,
                      format=None,
                      bands=None,
+                     combineBands=None,
                      srcNodata=None,
                      ovr=None,
                      targetCoordinateSystem=None,
@@ -3429,6 +3430,8 @@ def FootprintOptions(options=None,
         output format ("GeoJSON", etc...)
     bands:
         list of output bands to burn values into
+    combineBands:
+        how to combine bands: "union" (default) or "intersection"
     srcNodata:
         source nodata value(s).
     ovr:
@@ -3476,6 +3479,8 @@ def FootprintOptions(options=None,
         if bands is not None:
             for b in bands:
                 new_options += ['-b', str(b)]
+        if combineBands:
+            new_options += ["-combine_bands", combineBands]
         if targetCoordinateSystem:
             new_options += ["-t_cs", targetCoordinateSystem]
         if dstSRS: