Imitating MapBox' Cloudless Atlas

Inspired by several posts which explain MapBox' marvellous Cloudless Atlas, I wanted to imitate a cloudless atlas by myself. Limited in time and computer power I had to do it much simpler of course. But using first class spatial Python libraries I could achieve quite nice results with a simple Python script.

The result: a cloudless Landsat 8 scene

Armchair mapping using GRASS GIS

After reading my previous post a friend asked me why I didn't use just GRASS GIS. That's really a valid question. Since the GDAL/OGR utilities are some of my daily tools, I'm more familiar with them than with GRASS. Out of interest I developed a workflow how to extract and vectorize waterbodies from Landsat images using (almost solely) GRASS GIS. The resulting features I wanted to upload to OpenStreetMap.
This time the Nam Ngum 1 hydro power reservoir in Laos was my test area. Like Nam Leuk last time it was only roughly digitized and some of its islands needed to be remapped.

Garmin maps with contour lines

Some days ago I read in the OpenStreetMap wiki about maps for Garmin devices featuring contour lines. Interested in topographic Garmin maps for openstreetmap.la, I started to create maps with contour lines for Laos. As often happens it took me longer than I expected to figure out a suitable work flow, mainly because I ran in several memory issues.

I started with the generation of contour lines from the SRTM (and CGIAR improved) terrain model using gdal_contour. After about two hours calculating on a more capable work station (than my laptop), I got a 900MB heavy Shapefile. Based on this GDAL output, my plan was to convert the lines to the OSM format to use finally mkgmap.

To be able to deal reasonably with this amount of data I put the lines with shp2pgsql into a PostGIS database. Storing data in a PostGIS database facilitates faster spatial and attribute queries.

My first idea was to convert each contour level with ogr2osm to the OSM format. ogr2osm is a very flexible Python script I like to use, but in this case it refused to work for unknown reason. I had to use the two step approach with ogr2ogr and gpsbabel. For each contour level between 0 and 3000 I scripted the following steps and subdivided the levels in minor, medium and major contours according to the Garmin features land-contour-thin, land-contour-medium and land-contour-thick.

Query the database and write the result as tracks to a GPX file.

# ${level} means the currently processed contour level
ogr2ogr -f GPX -lco FORCE_GPX_TRACK=YES -sql "SELECT wkb_geometry, name FROM srtm_contours WHERE name=${level} AND ST_INTERSECTS(wkb_geometry, ST_SetSRID( ST_MakeBox2D( ST_POINT(100.1, 13.9),ST_POINT(107.7, 22.5)), 4326))" -nlt MULTILINESTRING ${outdir}/${prefix}_${level}.gpx PG:"dbname=openstreet user=openstreet"

Use gpsbabel to convert the GPX file to an OSM file and add corresponding tags.

gpsbabel -i gpx -f ${outdir}/${prefix}_${level}.gpx -o osm,tag="contour:elevation;contour_ext:${contour_ext}",created_by= -F ${outdir}/${prefix}_${level}_tmp.osm

Since it is not possible to add an elevation to a GPX track (only to points), I wrote the elevation to the name attribute. Now I had to replace the name tags with ele.

sed 's/name/ele/g'

The versatile osmosis program, I wanted to use to compress the data to PBF format, still didn't accept the file since it did not comply with the API v0.6. I had to add a version and timestamp attribute to each node and way and replace the negative identifiers with positive ones.

cat ${outdir}/${prefix}_${level}.osm | sed "s/node\ /node\ version='1'\ timestamp='2010-08-01T12:00:00+02:00' /g" | sed "s/way\ /way\ version='1'\ timestamp='2010-08-01T12:00:00+02:00' /g" | sed "s/'-/'/g" | osmosis --read-xml file=- --write-pbf omitmetadata=true ${outdir}/${prefix}_${level}.osm.pbf

Applying mkgmap to this file, mkgmap ran inevitable out of memory. It was necessary to split first the data (with splitter) to tiles manageable by mkgmap. I used the following tiles with dynamically generated map names for each tile and contour level:

KML split-file for splitter on Google Earth

Using this approach I finally managed to create map tiles in the Garmin IMG format.
Last step was to merge all contour lines maps with the OpenStreetMap data, but this could easily be done by mkgmap.

The result looks as follows:

Highest mountain in Laos

Cave near Vang Vieng

The map is can be downloaded as gmapsupp.img from openstreetmap.la or directly from here. There are detailed instructions on the OpenStreetMap wiki how to store a gmapsupp.img file in your Garmin device.

Update 9.9.2011: Due to problems with the routing, I've removed again the contours. To be continued ...