10

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. Large GTiff image

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:

enter image description here

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":

Image loaded into Cesium

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?

Improve this question
1
  • Why do you want to georeference the image? It already has all CRS information inside.
    – AndreJ
    Jun 28, 2016 at 15:51

2 Answers 2

Reset to default
7

I think all you need to do is reproject it using:

gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -t_srs EPSG:3857 image.tif newImage.tif

and then tile it:

gdal2tiles.py newImage.tif

If your file is very large it make take a while.

Improve this answer
8
  • You can control zoom levels with -z parameter, did you try it already? And notice that your image has 4 bands which may lead to unexpected results especially because you deal with 16 bit data. Some preprocessing may be needed first
    – user30184
    Jun 28, 2016 at 17:24
  • Thank you for your response @iant. I am going to try this and see what happens. Could you perhaps explain a bit more in your answer what the options you have selected mean? According to the documentation, -co option "passes a creation option to the output format driver". So are you in effect adding more properties to the tiff file?
    – effort
    Jun 29, 2016 at 6:54
  • @user30184 I did not try that yet no. How would I know what the zoom levels should be? Or can I specify them as I want? I thought that by leaving this option out, I would let the script determine the zoom level based on the area selected.
    – effort
    Jun 29, 2016 at 7:00
  • @user30184 you also mentioned that I might get unexpected results with a 4 banded image and 16 bit data. Why exactly? Is it not in the correct format for gdal tile processing? If so, the image was obtained directly from the provider, what would be the steps to get the file in the correct format? Namely a GDAL tiff file? (If I can say something like that)
    – effort
    Jun 29, 2016 at 7:02
  • I got the following problems running gdalwarp as you suggested @iant: >>>>>>>>>>>> ERROR 1: TIFFFillTile:Read error at row 43520, col 47104; got 35788250 bytes, expected 37421449 ERROR 1: TIFFReadEncodedTile() failed. ERROR 1: pleiades_merge05m_2015-06-19.tif, band 1: IReadBlock failed at X offset 86, Y offset 109 ERROR 1: GetBlockRef failed at X block offset 86, Y block offset 109 >>>>>>>>>>>> Any suggestions to fix this?
    – effort
    Jun 29, 2016 at 7:45
2

I suppose that your image is some of the the 4-band products by Airbus DS:

http://www.intelligence-airbusds.com/en/4951-which-spectral-mode-do-i-choose

Gdal2tiles is made for splitting common, visual images into png tiles. Such images use 8 bits per band and they have one band (greyscale), 3 bands (red-green-blue) of 4 bands (reg-green-blue + alpha).

I would say that your question is to large extent irrelevant because your source data do not suit well for gdal2tiles. You may get around the immediate problems you have right now but still the final result will still not be good if you do not re-process your data.

The reason for not-good-looking tile that you attached to your question may be that the fourth data band is interpreted to be an alpha channel.

Improve this answer
3
  • Thank you @user30184, I have been reading up on some resources and have come to a similar thought. I think it would be best to ask our providers to provide us with "GDAL compatible" tiff files, but until they get back to us, this is all we have. I was considering perhaps removing one of the bands to avoid gdal to interperting the last band as an alpha channel. Any suggestions?
    – effort
    Jun 29, 2016 at 11:44
  • Use gdal_translate for cutting a small subset from the image gdal_translate -srcwin 20000 20000 1000 1000 original.tif sample.tifshould do it (large offsets for avoiding the nodata areas). Open this small image with QGIS and you should be able to play with the viewing settings fast. My guess about the alpha channel is probably wrong, otherwise the result should look colorful, not grey.
    – user30184
    Jun 29, 2016 at 12:02
  • Thanks @user30184, I did as you suggested and successfully opened it in QGIS. See this link: drive.google.com/open?id=0B97NtaPJrVz-anRYQmxjZFludk0 How can I now go about "debugging" my issue? Using QGIS to do the gdalwarp and gdal2tiles?
    – effort
    Jun 29, 2016 at 13:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.