gdal2tiles: add --special-values and --special-value-pct-threshold sw…

…itches

autotest/pyscripts/test_gdal2tiles.py

@@ -33,12 +33,13 @@
 import os
 import os.path
 import shutil
+import struct
 import sys
 
 import pytest
 import test_py_scripts  # noqa  # pylint: disable=E0401
 
-from osgeo import gdal  # noqa
+from osgeo import gdal, osr  # noqa
 from osgeo_utils.gdalcompare import compare_db
 
 pytestmark = pytest.mark.skipif(
@@ -581,3 +582,44 @@ def test_gdal2tiles_py_webp(script_path, tmp_path, resampling):
         )
         diff_found = compare_db(gdal.Open(webp_filename), gdal.Open(filename))
         assert not diff_found, (resampling, filename)
+
+
+@pytest.mark.require_driver("PNG")
+def test_gdal2tiles_special_values(script_path, tmp_path):
+
+    input_tif = str(tmp_path / "test_gdal2tiles_special_values.tif")
+    output_folder = str(tmp_path / "test_gdal2tiles_special_values")
+
+    src_ds = gdal.GetDriverByName("GTiff").Create(input_tif, 256, 256, 3, gdal.GDT_Byte)
+    src_ds.GetRasterBand(1).WriteRaster(
+        0, 0, 2, 2, struct.pack("B" * 4, 10, 20, 30, 40)
+    )
+    src_ds.GetRasterBand(2).WriteRaster(
+        0, 0, 2, 2, struct.pack("B" * 4, 11, 21, 31, 41)
+    )
+    src_ds.GetRasterBand(3).WriteRaster(
+        0, 0, 2, 2, struct.pack("B" * 4, 12, 22, 32, 42)
+    )
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(3857)
+    src_ds.SetSpatialRef(srs)
+    MAX_GM = 20037508.342789244
+    RES_Z0 = 2 * MAX_GM / 256
+    RES_Z1 = RES_Z0 / 2
+    # Spatial extent of tile (0,0) at zoom level 1
+    src_ds.SetGeoTransform([-MAX_GM, RES_Z1, 0, MAX_GM, 0, -RES_Z1])
+    src_ds = None
+
+    test_py_scripts.run_py_script_as_external_script(
+        script_path,
+        "gdal2tiles",
+        f"-q -z 0-1 --special-values=30,31,32 --special-value-pct-threshold=50 {input_tif} {output_folder}",
+    )
+
+    ds = gdal.Open(f"{output_folder}/0/0/0.png")
+    assert struct.unpack("B" * 4, ds.ReadRaster(0, 0, 1, 1)) == (
+        (10 + 20 + 40) // 3,
+        (11 + 21 + 41) // 3,
+        (12 + 22 + 42) // 3,
+        255,
+    )

doc/source/programs/gdal2tiles.rst

@@ -22,6 +22,7 @@ Synopsis
                   [-w <webviewer>] [-t <title>] [-c <copyright>]
                   [--processes=<NB_PROCESSES>] [--mpi] [--xyz]
                   [--tilesize=<PIXELS>] [--tmscompatible]
+                  [--special-values=<SPECIAL_VALUES>] [--special-value-pct-threshold=<SPECIAL_VALUE_PCT_THRESHOLD>]
                   [-g <googlekey] [-b <bingkey>] <input_file> [<output_dir>] [<COMMON_OPTIONS>]
 
 Description
@@ -143,6 +144,23 @@ can publish a picture without proper georeferencing too.
 
   .. versionadded:: 3.6
 
+.. option:: --special-values=<SPECIAL_VALUES>
+
+  Comma-separated tuple of values (thus typically "R,G,B"), that are ignored
+  as contributing source * pixels during resampling. The number of values in
+  the tuple must be the same as the number of bands, excluding the alpha band.
+  Several tuples of special values may be specified using the "(R1,G1,B2),(R2,G2,B2)" syntax.
+  Only taken into account by Average currently.
+  This concept is a bit similar to nodata/alpha, but the main difference is
+  that pixels matching one of the special value tuples are still considered
+  as valid, when determining the target pixel validity/density.
+
+.. option:: --special-value-pct-threshold=SPECIAL_VALUE_PCT_THRESHOLD
+
+  Minimum percentage of source pixels that must be set at one of the --special-values to cause the special
+  value, that is in majority among source pixels, to be used as the target pixel value. Default value is 50(%)
+
+  .. versionadded:: 3.9
 
 .. option:: -h, --help
 

swig/python/gdal-utils/osgeo_utils/gdal2tiles.py

@@ -873,7 +873,27 @@ def scale_query_to_tile(dsquery, dstile, options, tilefilename=""):
     tile_size = dstile.RasterXSize
     tilebands = dstile.RasterCount
 
-    if options.resampling == "average":
+    dsquery.SetGeoTransform(
+        (
+            0.0,
+            tile_size / float(querysize),
+            0.0,
+            0.0,
+            0.0,
+            tile_size / float(querysize),
+        )
+    )
+    dstile.SetGeoTransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0))
+
+    if options.resampling == "average" and options.special_values:
+
+        gdal.Warp(
+            dstile,
+            dsquery,
+            options=f"-r average -wo SPECIAL_VALUES={options.special_values} -wo SPECIAL_VALUE_PCT_THRESHOLD={options.special_value_pct_threshold}",
+        )
+
+    elif options.resampling == "average":
 
         # Function: gdal.RegenerateOverview()
         for i in range(1, tilebands + 1):
@@ -949,18 +969,6 @@ def scale_query_to_tile(dsquery, dstile, options, tilefilename=""):
             gdal_resampling = gdal.GRA_Q3
 
         # Other algorithms are implemented by gdal.ReprojectImage().
-        dsquery.SetGeoTransform(
-            (
-                0.0,
-                tile_size / float(querysize),
-                0.0,
-                0.0,
-                0.0,
-                tile_size / float(querysize),
-            )
-        )
-        dstile.SetGeoTransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0))
-
         res = gdal.ReprojectImage(dsquery, dstile, None, None, gdal_resampling)
         if res != 0:
             exit_with_error(
@@ -1346,6 +1354,7 @@ def create_base_tile(tile_job_info: "TileJobInfo", tile_detail: "TileDetail") ->
     # Tile dataset in memory
     tilefilename = os.path.join(output, str(tz), str(tx), "%s.%s" % (ty, tileext))
     dstile = mem_drv.Create("", tile_size, tile_size, tilebands)
+    dstile.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand)
 
     data = alpha = None
 
@@ -1399,6 +1408,8 @@ def create_base_tile(tile_job_info: "TileJobInfo", tile_detail: "TileDetail") ->
             # Big ReadRaster query in memory scaled to the tile_size - all but 'near'
             # algo
             dsquery = mem_drv.Create("", querysize, querysize, tilebands)
+            dsquery.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand)
+
             # TODO: fill the null value in case a tile without alpha is produced (now
             # only png tiles are supported)
             dsquery.WriteRaster(
@@ -1422,6 +1433,11 @@ def create_base_tile(tile_job_info: "TileJobInfo", tile_detail: "TileDetail") ->
             tilefilename, dstile, strict=0, options=_get_creation_options(options)
         )
 
+        # Remove useless side car file
+        aux_xml = tilefilename + ".aux.xml"
+        if gdal.VSIStatL(aux_xml) is not None:
+            gdal.Unlink(aux_xml)
+
     del dstile
 
     # Create a KML file for this tile.
@@ -1485,10 +1501,12 @@ def create_overview_tile(
     dsquery = mem_driver.Create(
         "", 2 * tile_job_info.tile_size, 2 * tile_job_info.tile_size, tilebands
     )
+    dsquery.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand)
     # TODO: fill the null value
     dstile = mem_driver.Create(
         "", tile_job_info.tile_size, tile_job_info.tile_size, tilebands
     )
+    dstile.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand)
 
     usable_base_tiles = []
 
@@ -1759,6 +1777,19 @@ def optparse_init() -> optparse.OptionParser:
         type="choice",
         help="which tile driver to use for the tiles",
     )
+    p.add_option(
+        "--special-values",
+        dest="special_values",
+        type=str,
+        help="Tuples of values (e.g. <R>,<G>,<B> or (<R1>,<G1>,<B1>),(<R2>,<G2>,<B2>)) that must be ignored as contributing source pixels during resampling. Only taken into account for average resampling",
+    )
+    p.add_option(
+        "--special-value-pct-threshold",
+        dest="special_value_pct_threshold",
+        type=float,
+        default=50,
+        help="Minimum percentage of source pixels that must be set at one of the --special-values to cause the special value, that is in majority among source pixels, to be used as the target pixel value. Default value is 50 (%)",
+    )
 
     # KML options
     g = optparse.OptionGroup(