Blue Marble Mosaic

Recently, I decided to push www.olympianengine.com to Google App Engine. Google App Engine is a platform for building and hosting web applications on Google Web Servers (wikipedia).

I hit an immediate show-stopper: the maximum size for any single file is 1048576 bytes. I spent a night chopping up the bluemarble.tif from the GeoTIFF tutorial, such that I can load it in chunks, then realized I have large KML files on another website that I have no intention of chopping up.

Thus, Google App Engine, does not have the capacity to support my KML files, and is no longer a viable solution for me at all. However, out of the ashes, comes this tutorial on how to use the Geospatial Data Abstraction Library (GDAL) utilities to manipulate GeoTIFF images.

Goals:
a. Build GDAL on local system.
b. Use gdalwarp to manually cut up a large GeoTIFF into smaller pieces.
c. Use gdal_merge.py to put together the set of files we made in (b).

Environment:
Ubuntu 8.10
GNU Make version 3.81
gcc version 4.3.2

1/4 Build GDAL on local system.
Download GDAL:
ftp://ftp.remotesensing.org/gdal/

From /home/myuser/:

$ gzip -rd gdal-1.5.2.tar.gz
$ tar -xf gdal-1.5.2.tar
$ cd gdal-1.5.2
$ mkdir deploy
$ ./configure --prefix=/home/myuser/gdal-1.5.2/deploy --with-python
$ make
$ make install

Additional information for building GDAL with Python:
http://trac.osgeo.org/gdal/wiki/BuildHints
http://trac.osgeo.org/gdal/wiki/GdalOgrInPython

Make sure to change the following environment variables after a successful build:

export PATH=$PATH:/home/myuser/gdal-1.5.2/deploy/bin
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/myuser/gdal-1.5.2/deploy/lib
export PYTHONPATH=$PYTHONPATH:/home/myuser/gdal-1.5.2/deploy/lib/python2.5/site-packages

Now, we’ll be able to use the GDAL utilities from the command line.

2/4 Use gdalwarp to cut bluemarble.tif:
Download bluemarble.tif.

We’re going to use gdalwarp to cut up the bluemarble.tif into 32 seperate files, our tile grid can be seen in the below image.

bluemarblegrid

Documentation for gdalwarp can be seen here:
http://www.gdal.org/gdalwarp.html

Our command will take the form of: $ gdalwarp -te xmin ymin xmax ymax inputfile outputfile

$ gdalwarp -te -180 45 -135 90 bluemarble.tif 0.tif
$ gdalwarp -te -135 45 -90 90 bluemarble.tif 1.tif
$ gdalwarp -te -90 45 -45 90 bluemarble.tif 2.tif
$ gdalwarp -te -45 45 0 90 bluemarble.tif 3.tif
$ gdalwarp -te 0 45 45 90 bluemarble.tif 4.tif
$ gdalwarp -te 45 45 90 90 bluemarble.tif 5.tif
$ gdalwarp -te 90 45 135 90 bluemarble.tif 6.tif
$ gdalwarp -te 135 45 180 90 bluemarble.tif 7tif
$ gdalwarp -te -180 0 -135 45 bluemarble.tif 8.tif
$ gdalwarp -te -135 0 -90 45 bluemarble.tif 9.tif
$ gdalwarp -te -90 0 -45 45 bluemarble.tif 10.tif
$ gdalwarp -te -45 0 0 45 bluemarble.tif 11.tif
$ gdalwarp -te 0 0 45 45 bluemarble.tif 12.tif
$ gdalwarp -te 45 0 90 45 bluemarble.tif 13.tif
$ gdalwarp -te 90 0 135 45 bluemarble.tif 14.tif
$ gdalwarp -te 135 0 180 45 bluemarble.tif 15.tif
$ gdalwarp -te -180 -45 -135 0 bluemarble.tif 16.tif
$ gdalwarp -te -135 -45 -90 0 bluemarble.tif 17.tif
$ gdalwarp -te -90 -45 -45 0 bluemarble.tif 18.tif
$ gdalwarp -te -45 -45 0 0 bluemarble.tif 19.tif
$ gdalwarp -te 0 -45 45 0 bluemarble.tif 20.tif
$ gdalwarp -te 45 -45 90 0 bluemarble.tif 21.tif
$ gdalwarp -te 90 -45 135 0 bluemarble.tif 22.tif
$ gdalwarp -te 135 -45 180 0 bluemarble.tif 23.tif
$ gdalwarp -te -180 -90 -135 -45 bluemarble.tif 24.tif
$ gdalwarp -te -135 -90 -90 -45 bluemarble.tif 25.tif
$ gdalwarp -te -90 -90 -45 -45 bluemarble.tif 26.tif
$ gdalwarp -te -45 -90 0 -45 bluemarble.tif 27.tif
$ gdalwarp -te 0 -90 45 -45 bluemarble.tif 28.tif
$ gdalwarp -te 45 -90 90 -45 bluemarble.tif 29.tif
$ gdalwarp -te 90 -90 135 -45 bluemarble.tif 30.tif
$ gdalwarp -te 135 -90 180 -45 bluemarble.tif 31.tif

We now have 32 tiles, each underneath the target size of 1048576 bytes. Each tile is also properly geospatially enabled, which we can see when we run gdalinfo on the first image. If someone was so inclined they could make a simple shell script with two loops that would effectively do all the above in a few lines.

$ gdalinfo 0.tif
Driver: GTiff/GeoTIFF
Files: 0.tif
Size is 540, 540
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235630016,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.083333333333332,-0.083333333333332)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0'0.00"W, 90d 0'0.00"N)
Lower Left  (-180.0000000,  45.0000000) (180d 0'0.00"W, 45d 0'0.00"N)
Upper Right (-135.0000000,  90.0000000) (135d 0'0.00"W, 90d 0'0.00"N)
Lower Right (-135.0000000,  45.0000000) (135d 0'0.00"W, 45d 0'0.00"N)
Center      (-157.5000000,  67.5000000) (157d30'0.00"W, 67d30'0.00"N)
Band 1 Block=540x5 Type=Byte, ColorInterp=Red
Band 2 Block=540x5 Type=Byte, ColorInterp=Green
Band 3 Block=540x5 Type=Byte, ColorInterp=Blue

Coordinate system, coordinate corners as well as pixel size indicate that each tile was created properly.

Files: gdalwarpedfiles.tar.gz

3/4 Use gdal_merge.py to put the tiles back into one cohesive bluemarble.tif:

$ gdal_merge.py -o ../bluemarblemerged.tif *.tif

We now have a complete reconstruction of our original GeoTIFF. Running gdalinfo on the the image shows us the geospatial information commiserate with the original bluemarble.tif.

$ gdalinfo ../bluemarblemerge.tif
Driver: GTiff/GeoTIFF
Files: ../bluemarblemerge.tif
Size is 4320, 2160
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235630016,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.083333333333332,-0.083333333333332)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0'0.00"W, 90d 0'0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0'0.00"W, 90d 0'0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0'0.00"E, 90d 0'0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0'0.00"E, 90d 0'0.00"S)
Center      (  -0.0000000,   0.0000000) (  0d 0'0.00"W,  0d 0'0.00"N)
Band 1 Block=4320x1 Type=Byte, ColorInterp=Red
Band 2 Block=4320x1 Type=Byte, ColorInterp=Green
Band 3 Block=4320x1 Type=Byte, ColorInterp=Blue

4/4 Complete:
Possible extensions:
a) Creating a script to run the gdalwarp commands based on how many tiles in each x/y direction are desired.
b) Use gdal2tiles.py to create a tile set.

Leave a Reply