From 861d286b8c7656f16715017cea39e0368d1b02a1 Mon Sep 17 00:00:00 2001 From: Even Rouault 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 3f5e21118630..d7e996c0be0b 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 2a463cc8b7eaa65d44d28b9b1c8d7f508f3a202a Mon Sep 17 00:00:00 2001 From: Even Rouault 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 3c516452c615cb010d5105ca0a6c66d97f150084 Mon Sep 17 00:00:00 2001 From: Even Rouault 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 d7e996c0be0b..68f86f1d0d1a 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)))")