[GIS] How to use gdal2tiles on a custom tiff image recieved from providers to generate tiles

cesiumgdalgdal2tilesgdalwarp

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?

Best Answer

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.

Related Solutions

[GIS] Georeferencing png using gdal

Here's a world file for you (as per Gene's suggestion). Cut and paste the following into a text editor and name the file the same as your full-sized map but give it the extension '.pgw' instead on '.png':

1.0
0.0
0.0
-1.0
0.0
3231.0

I've made the following assumptions:

Your map has a resolution (pixel size) of 1m square (i.e. each pixel is a reasonable stride length for a fantasy game humanoid character). If I am wrong, edit lines 1 and 4 accordingly.
I have assumed that your map is not rotated.
I have assumed that your LL corner is at 0,0 in your world. If your pixel size is bigger than 1m square then you will need to multiply your Y coordinate (3231.0) accordingly.
Technically the UL coordinate (last two lines) is the center of your pixel. I have been lazy and set it to the corner. This in theory means that your coordinates would be consistently half a meter out but it's a fantasy world so nobody will know; it is consistent so won't screw up any distant calculations; and it saves a negative value as one of your origin coordinates. Adjust if you are sufficiently worried :)

I would then re-tile your map from the georeferenced full-sized one (this saves you making a world file for each tile)

Now, if you load your fantasy map into a GIS or similar system, you will probably be asked for the Spatial Reference System (SRS). The World file doesn't encapsulate that. It just locates and orients the image file in space relative to a nominal 0,0 origin. Since this is a fantasy map, all you need to do is pretend that it is in the real world somewhere and pick any coordinate system that gives you a flat projection with units in meters (or feet if you prefer but you will may have to adjust the values for resolution in the world file). A couple of suggestions for SRS could be:

EPSG: 3857 (aka 900913)
EPSG: 27700 (UK National Grid - the advantage of a national grid system UK or otherwise, is that you mostly stay clear of the poles and the international dateline and can hopefully avoid any weirdness which may result from such a location)

Your choice might not matter at all... on the other hand, you choice of SRS might be influenced by what you intend to do with your KML/CSV file (like you might be making a game with a post-apocalyptic theme but set in Quebec ...with zombies...). If you are using Google Earth or similar as a 'platform' for your game/fantasy world, so long as all your other data are in the same coordinate system, all should be well.

However, I do wonder if you may be better sticking with your PNG file as a text format like KML/CSV will be big. A PNG may be better for an internet-based game (reducing bandwidth). If you have a stand-alone scenario, most programming languages understand image formats well (especially if you call the GDAL libraries) and so I'm not convinced that a text-based format is even necessary here - but that's your choice of course.

[GIS] Can gdalwarp reproject from ESPG:4326 (WGS84) to UTM

Try with other zone it has no problem, I guess GDAL cannot wrap it outside.

gdalwarp -t_srs '+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs' test.tif test_utm1.tif

Best Answer

Related Solutions

[GIS] Georeferencing png using gdal

[GIS] Can gdalwarp reproject from ESPG:4326 (WGS84) to UTM

Related Question