From 6fd91c848dfdb879a8f736d51ed65cf4692500ac Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 27 Nov 2023 15:33:13 +0100
Subject: [PATCH 1/3] gdal_footprint: really fix datasets with alpha band, by
 vectorizing only the alpha band and not the other ones (really fixes #8792,
 fixes #8834)

---
 apps/gdal_footprint_lib.cpp                   | 3 +--
 autotest/utilities/test_gdal_footprint_lib.py | 3 +--
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/apps/gdal_footprint_lib.cpp b/apps/gdal_footprint_lib.cpp
index aa09b6d378a8..ae93a7ade69a 100644
--- a/apps/gdal_footprint_lib.cpp
+++ b/apps/gdal_footprint_lib.cpp
@@ -703,8 +703,7 @@ 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;
             }
diff --git a/autotest/utilities/test_gdal_footprint_lib.py b/autotest/utilities/test_gdal_footprint_lib.py
index 8accaa983330..999770b67438 100755
--- 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,

From c95f6e725341f3e13575904fd6714a6e4f987302 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 27 Nov 2023 16:08:25 +0100
Subject: [PATCH 2/3] Python bindings: add a combineBands option to
 gdal.Footprint()

---
 swig/include/python/gdal_python.i | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/swig/include/python/gdal_python.i b/swig/include/python/gdal_python.i
index d49c2130b0d9..54becf4b2f4e 100644
--- 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:

From 15e3168a226a8cb8818a190f0cda097ab2858c32 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 27 Nov 2023 16:08:57 +0100
Subject: [PATCH 3/3] gdal_footprint: fix taking into account of individual
 bands that have nodata

---
 apps/gdal_footprint_lib.cpp                   |  2 +-
 autotest/utilities/test_gdal_footprint_lib.py | 83 +++++++++++++++++++
 2 files changed, 84 insertions(+), 1 deletion(-)

diff --git a/apps/gdal_footprint_lib.cpp b/apps/gdal_footprint_lib.cpp
index ae93a7ade69a..10e88c224c4a 100644
--- a/apps/gdal_footprint_lib.cpp
+++ b/apps/gdal_footprint_lib.cpp
@@ -709,7 +709,7 @@ static bool GDALFootprintProcess(GDALDataset *poSrcDS, OGRLayer *poDstLayer,
             }
             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
index 999770b67438..84b8760d6e89 100755
--- a/autotest/utilities/test_gdal_footprint_lib.py
+++ b/autotest/utilities/test_gdal_footprint_lib.py
@@ -451,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)))")