Fixed gdal2tiles.py alpha channel handling

[trinopoty]

What does this PR do?

gdal2tiles.py by default prefers to use mask band over alpha band.
In certain edge cases, the mask indicates all data as valid and the alpha band contains the actual transparency information. This results in gdal2tiles.py creating output files with black regions.
This PR attempts to fix this case and force the use of an alpha band if it's present, otherwise, fall back to a mask band.

What are related issues/pull requests?

Tasklist

• Add test case(s)
• Review
• Adjust for comments
• All CI builds and checks have passed

Environment

Unrelated

Any advice on how to implement tests would be appreciated.

[trinopoty]

I encountered this problem while trying to process output from Pix4D.

[rouault]

can you show the output of gdalinfo on input files that would require that patch ?

[trinopoty]

Here is one such file:

Driver: GTiff/GeoTIFF
Files: ortho.tif
Size is 63946, 80189
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 44N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 44N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",81,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - 78°E to 84°E - by country"],
        BBOX[0,78,84,84]],
    ID["EPSG",32644]]
Data axis to CRS axis mapping: 1,2
Origin = (624084.719980000052601,2406254.382970000151545)
Pixel Size = (0.028280000000000,-0.028280000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  624084.720, 2406254.383) ( 82d12' 0.07"E, 21d45'20.06"N)
Lower Left  (  624084.720, 2403986.638) ( 82d11'59.46"E, 21d44' 6.32"N)
Upper Right (  625893.113, 2406254.383) ( 82d13' 3.02"E, 21d45'19.60"N)
Lower Right (  625893.113, 2403986.638) ( 82d13' 2.40"E, 21d44' 5.86"N)
Center      (  624988.916, 2405120.511) ( 82d12'31.23"E, 21d44'42.96"N)
Band 1 Block=63946x1 Type=Byte, ColorInterp=Red
  NoData Value=-10000
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=63946x1 Type=Byte, ColorInterp=Green
  NoData Value=-10000
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=63946x1 Type=Byte, ColorInterp=Blue
  NoData Value=-10000
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=63946x1 Type=Byte, ColorInterp=Alpha
  NoData Value=-10000

[rouault]

I'm a bit skeptical about your change. On the file you exhibit, GetMaskBand() and GetBand(alpha_band_number) should both return Band 4

[trinopoty]

I could not check if GetMaskBand() is returning Band 4 because the python API does not have a method or property for getting the band number from a Band object that I'm aware of.
The unpatched version of gdal2tiles.py creates loads of black filled png files and the patch fixes that.

[trinopoty]

Let me try to dump out the mask and the alpha band into separate files.

[rouault]

the python API does not have a method or property for getting the band number from a Band object that I'm aware of.

==> band.GetBand() (maps to GDALGetBandNumber())

[trinopoty]

Both GetMaskBand() and GetBand(alpha_band_number) does seem to refer to the 4th band.
But then, why does GetMaskBand() produce black images? And how does this patch fix that if both are equivalent?

[rouault]

And how does this patch fix that if both are equivalent?

that's why I can't make sense of it. Are you sure you get different results with and without your patch on the very image you mentioned in #2299 (comment) and not on another one ?

[trinopoty]

You can see the two different output produced by gdal2tiles.py.
One with the patch and one without the patch on the exact same input file.

[trinopoty]

If you want, I could try to upload the file somewhere and link it but it's 6GB in size.

[rouault]

Wait, there's one subtlety I overlooked. Your file has both nodata set (to -10000) and an alpha band (band 4). The GetMaskBand() logic would normally return a synthetic mask band based on the nodata value, but as -10000 is out of range for Byte data type, it goes on to select an alpha band. What you observed. But perhaps you've not instrumented all points, and there's some transformations within gdal2tiles that change the settings a bit.

I could try to upload the file somewhere and link it but it's 6GB in size.

Try to minimize it first with gdal_translate -outsize 1% 1% for example and see if you can still reproduce the issue

[trinopoty]

I was able to reproduce the bug with the resized tif.
Here is the tif file with the output from the patched and unpatched gdal2tiles.py.

[trinopoty]

But perhaps you've not instrumented all points, and there's some transformations within gdal2tiles that change the settings a bit.

I'm not all that familiar with the internals of gdal and it's entirely possible I'm overlooking something.

[rouault]

gdal2tiles is confused by the invalid range of the nodata value. Try the following patch that will detect that

diff --git a/gdal/swig/python/scripts/gdal2tiles.py b/gdal/swig/python/scripts/gdal2tiles.py
index 65288d1..5e2ada4 100755
--- a/gdal/swig/python/scripts/gdal2tiles.py
+++ b/gdal/swig/python/scripts/gdal2tiles.py
@@ -695,8 +695,13 @@ def setup_no_data_values(input_dataset, options):
             in_nodata = nds
     else:
         for i in range(1, input_dataset.RasterCount + 1):
-            raster_no_data = input_dataset.GetRasterBand(i).GetNoDataValue()
+            band = input_dataset.GetRasterBand(i)
+            raster_no_data = band.GetNoDataValue()
             if raster_no_data is not None:
+                if band.DataType == gdal.GDT_Byte and (raster_no_data != int(raster_no_data) or raster_no_data < 0 or raster_no_data > 255):
+                    # We should possibly do similar check for other data types
+                    in_nodata = []
+                    break
                 in_nodata.append(raster_no_data)
 
     if options.verbose:

I'm wondering if a better/complementary fix wouldn't be for the GTiff driver to not report such non sensical nodata value to begin with (with a warning)

[rouault]

I've committed the above patch in 9a16f92. Should fix the issue

[trinopoty]

I have tested and this issue has been resolved by your commit (34f25f0).