Quality of gdal2tiles output

Question

I am trying to generate a tileset from a GeoTIFF using gdal2tiles, e.g.:

gdal2tiles -s EPSG:3448 input.tif outputfolder

But the output tiles’ quality is definitely not what I expected (see images below).

Looking at the docs I don’t see which parameter could fix the issue. I played with the resampling methods but the result does not look any better.

Any idea how to go about this or what other tools to use before generating the tiles?

Input:

Output:

Input:

Output:

EDIT:

Here's the result of gdalinfo (the image is called 11e.tif):

Driver: GTiff/GeoTIFF
Files: 11e.tif
       11e.tif.ovr
       11e.tif.aux.xml
Size is 20000, 20000
Coordinate System is:
PROJCS["JAD_2001_Jamaica_Grid",
    GEOGCS["GCS_JAD_2001",
        DATUM["Jamaica_2001",
            SPHEROID["WGS_1984",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6758"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",18],
    PARAMETER["central_meridian",-77],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",750000],
    PARAMETER["false_northing",650000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (760000.000000000000000,650000.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  760000.000,  650000.000) ( 76d54'20.07"W, 17d59'59.92"N)
Lower Left  (  760000.000,  640000.000) ( 76d54'20.25"W, 17d54'34.65"N)
Upper Right (  770000.000,  650000.000) ( 76d48'40.15"W, 17d59'59.67"N)
Lower Right (  770000.000,  640000.000) ( 76d48'40.49"W, 17d54'34.41"N)
Center      (  765000.000,  645000.000) ( 76d51'30.24"W, 17d57'17.18"N)
Band 1 Block=128x128 Type=UInt16, ColorInterp=Gray
  Min=0.000 Max=2047.000 
  Minimum=0.000, Maximum=2047.000, Mean=829.688, StdDev=247.897
  NoData Value=65536
  Overviews: 10000x10000, 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313, 157x157
  Metadata:
    STATISTICS_COVARIANCES=61452.9798540855,71240.91544194205,80354.57525810145,50760.91993708123
    STATISTICS_MAXIMUM=2047
    STATISTICS_MEAN=829.688137455
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=247.89711546141
Band 2 Block=128x128 Type=UInt16, ColorInterp=Undefined
  Min=0.000 Max=2047.000 
  Minimum=0.000, Maximum=2047.000, Mean=818.696, StdDev=301.653
  NoData Value=65536
  Overviews: 10000x10000, 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313, 157x157
  Metadata:
    STATISTICS_COVARIANCES=71240.91544194205,90994.83133975854,102183.725048214,77838.12949558289
    STATISTICS_MAXIMUM=2047
    STATISTICS_MEAN=818.6958275225
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=301.65349548739
Band 3 Block=128x128 Type=UInt16, ColorInterp=Undefined
  Min=0.000 Max=2047.000 
  Minimum=0.000, Maximum=2047.000, Mean=545.815, StdDev=350.652
  NoData Value=65536
  Overviews: 10000x10000, 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313, 157x157
  Metadata:
    STATISTICS_COVARIANCES=80354.57525810145,102183.725048214,122957.1451096333,93612.24267034492
    STATISTICS_MAXIMUM=2047
    STATISTICS_MEAN=545.8145240525
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=350.65245630058
Band 4 Block=128x128 Type=UInt16, ColorInterp=Undefined
  Min=0.000 Max=2047.000 
  Minimum=0.000, Maximum=2047.000, Mean=388.913, StdDev=400.207
  NoData Value=65536
  Overviews: 10000x10000, 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313, 157x157
  Metadata:
    STATISTICS_COVARIANCES=50760.91993708123,77838.12949558289,93612.24267034492,160165.495608783
    STATISTICS_MAXIMUM=2047
    STATISTICS_MEAN=388.91275769
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=400.20681604488

EDIT 2:

According to the gdal2tiles docs:

Inputs with non-Byte data type (i.e. Int16, UInt16,…) will be clamped to the Byte data type, causing wrong results. To awoid this it is necessary to rescale input to the Byte data type using gdal_translate utility.

But if I run the following I get a strange looking dataset:

gdal_translate 11e.tif -ot Byte -of GTIFF -b 1 -b 2 -b 3 11e_byte_multiband.tif

Needless to say the tiles generated from this output look just as strange.

With a couple more parameters it looks better but still off from the original:

gdal_translate -ot Byte -of GTIFF -scale -co tiled=yes -co compress=deflate -b 1 -b 2 -b 3 11e.tif 11e_byte_multiband.tif

Answer 1

As you noticed from the gdal2tiles documentation https://gdal.org/programs/gdal2tiles.html, the Python script deals well only with 8-bit input images.

Note

Inputs with non-Byte data type (i.e. Int16, UInt16,…) will be clamped to the Byte data type, causing wrong results. To awoid this it is necessary to rescale input to the Byte data type using gdal_translate utility.

In addition your input imagery is 4-band while gdal2tiles.py is made for 1 or 3 data bands. Therefore you must use gdal_translate, select the three bands to use for RGB, and scale the band values to fit nicely into 8 bit. The command to use is more or less like

gdal_translate -ot Byte -of GTIFF -scale -co tiled=yes -co compress=deflate -b 1 -b 2 -b 3 input_16_bit_multiband.tif output_8_bit_RGB.tif

This command creates also a tiled GeoTIFF wiht lossless compression that yields a smaller and better manageable file for furthrer processing. If the histogram of the source image is very special you may need to give the input and output ranges for -scale manually. See the syntax from https://gdal.org/programs/gdal_translate.html.