I have been struggling a bit to generate tiles for a high-resolution image that we have. The current image we have is a very large (+20GB) image, saved as a GeoTiff file.
I would like to generate the tiles using the gdal2tiles command line utility and then open and view it in Cesium, using the TMS imagery provider to provide the tiles. Using gdalinfo, here is some of the details of the image:
Driver: GTiff/GeoTIFF
Files: image.tif
Size is 52250, 56119
Coordinate System is:
PROJCS["WGS 84 / UTM zone 35S",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",27],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",10000000],
UNIT["meters",1],
AUTHORITY["EPSG","32735"]]
Origin = (606276.000000000000000,7197873.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_MAXSAMPLEVALUE=13165
TIFFTAG_MINSAMPLEVALUE=1
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_SOFTWARE=ERDAS IMAGINE
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 606276.000, 7197873.000) ( 28d 3'21.59"E, 25d19'55.12"S)
Lower Left ( 606276.000, 7169813.500) ( 28d 3'29.55"E, 25d35' 7.17"S)
Upper Right ( 632401.000, 7197873.000) ( 28d18'55.92"E, 25d19'47.60"S)
Lower Right ( 632401.000, 7169813.500) ( 28d19' 5.85"E, 25d34'59.57"S)
Center ( 619338.500, 7183843.250) ( 28d11'13.23"E, 25d27'27.58"S)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray
Band 2 Block=512x512 Type=UInt16, ColorInterp=Undefined
Band 3 Block=512x512 Type=UInt16, ColorInterp=Undefined
Band 4 Block=512x512 Type=UInt16, ColorInterp=Undefined
My first attempt was to use gdal_translate to georeference the image, and then use gdalwarp to change the projection to EPSG:3857, as required by Cesium (see the API reference)
gdal_translate -of VRT -a_srs EPSG:4326 -gcp 606275 7197875 28.055987 -25.331974 -gcp 606275 7169814 28.058200 -25.585326 -gcp 632400.5 7197875 28.31553 -25.329876 -gcp 632400.5 7169814 28.318286 -25.583209 image.tif newImage1.vrt
gdalwarp -of VRT -t_srs EPSG:3857 newImage1.vrt newImage2.vrt
However, I get many of the following errors:
ERROR 1: latitude or longitude exceeded limits
Another method I tried was to use gdal2tiles directly, and generating the tiles:
gdal2tiles.py image.tif
This created a folder with one subfolder (labeled 18) being the only zoom level at which tiles were created. However, the images that I get here are completely "wrong" and "blurry".
An example of one of the tiles:
Any suggestions for generating tiles for this image large image of a specific area using gdal2tiles so that I can load and view it in Cesium?
Update
So, after trying @iant's suggestion, I used the following commands:
gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=YES -t_srs EPSG:3857 image.tif newImage.tif
This command worked perfectly fine until the very end where I got the following error:
ERROR 1: TIFFFillTile:Read error at row 43520, col 47104; got 35788250 bytes, expected 37421449
Not sure what this error meant, I left it for the moment and still got a final image "newImage.tif", produced by the gdalwarp step. Using this I called gdal2tiles.py
gdal2tiles.py newImage.tif
This produced a folder with subfolders 10-18 (and not just one zoom level 18 as I got previously). It also reads perfectly fine into Cesium, without any console errors, but the image still looks "wrong":
I am considering my issue may be as @user30184 has suggested "...source data does not suit well for gdal2tiles." However until our provider is able to provide us something for use with gdal, this is all I have.
I was considering perhaps removing one of the bands to avoid gdal to interperting the last band as an alpha channel. Any suggestions?