For a little while now I’ve been manipulating the excellent USGS geotiffs to fed my PCT mapping fetish. Since I’ve worked out a good enough system I thought I might share it with the world.
First off you’ll need a bunch of stuff installed on your computer.
With those installed you’re ready to follow my process.
1. Remove the collar
The USGS maps that I’ve been using come with a ‘collar’, that bit around the outside of the actual map data that you’d have to cut off if you were printing them out and sticking them together.
Before cutting the collar off make sure to keep a copy of the original tiff because most image editing programs will lose the geo data when you press Save. In Photoshop I had to convert the image to RGB mode so that I could make the outside part transparent. It’s important not to change the size of the image. The embedded coordinates are for the corner of the tiff even though the map doesn’t reach there. Save the image and get back to the command line for step 2.
2. Restore the geospatial data
The new tiff you have just saved does not have the embedded coordinates that it started with. So long as you kept a copy of the original file you can reinstate the data using the following command and the gdalcopyproj program.
> gdalcopyproj.py original_file.tif new_file.tif
3. Create a Virtual Tile Set
If you’re combining multiple maps you’ll either need to merge them into one large file or create a virtual tile set like this
>gdalbuildvrt -srcnodata b4 -hidenodata merged.vrt *.tif
The “b4″ tells GDAL that band four of your images is the transparency layer, and the “hidenodata” tells it that when two maps overlap, hide the one that has no data (is transparent). Without that the transparent collars of the geotifs erase useful content from geotifs they join up with.
You can merge the files like this
> gdalwarp -co COMPRESS=LZW *.tif merged.tif
But it won’t do the compression until the end, so with USGS maps you’ll be generating files over 1GB even if you only have a few source maps. The virtual tile set method is definitely quicker and takes less disk space.
4. Cutting the tiles
The following command will take either your merged geotif file, or the virtual tile set as the first parameter and an output directory as the second (which will be created if it doesn’t exist).
>gdal2tiles.py merged.vrt tiles
The process can take quite a while, but it starts with the most detailed layer and each one after that is about 4 times faster than the previous one. When it’s complete you’ll have a full TMS tile stack ready for use on web map systems such as Google Maps and OpenLayers and my iOS app Tiled Maps. Next time I’ll post the scripts I use to optimise a tile set down from ~1.7Gb to 700Mb.