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 ...

Set up pgRouting with OpenStreetMap data

Another rainy season Saturday got me time for some GIS exercises. Inspired by Carson Farmer's excellent introduction to pgRouting I decided to set up a routing enabled PostGIS database with OpenStreetMap data and compare it with my GRASS GIS routing setup.
Contrary to Carson Farmer I didn't compile pgRouting from source but installed it from the ubuntugis-unstable repository

sudo apt-get install postgres-8.4-pgrouting postgres-8.4-pgrouting-dd postgres-8.4-pgrouting-tsp

and followed the instruction from the mentioned blog.

The osm2pgrouting script creates all necessary tables including table "public.ways", which is probably the most important one. The table schema without indexes and constraints looks like:

Table "public.ways"
    Column    |       Type       | Modifiers 
--------------+------------------+-----------
 gid          | integer          | not null
 class_id     | integer          | not null
 length       | double precision | 
 name         | character(200)   | 
 x1           | double precision | 
 y1           | double precision | 
 x2           | double precision | 
 y2           | double precision | 
 reverse_cost | double precision | 
 rule         | text             | 
 to_cost      | double precision | 
 osm_id       | integer          | 
 the_geom     | geometry         | 
 source       | integer          | 
 target       | integer          | 

The relevant columns for the routing are "length","reverse_cost" and "to_cost". The latter two are the forward and backward costs to use a way and are used in directed shortest path calculations. It's necessary to use an algorithm that can handle direction when one-ways have to be taken in account.

Since OpenStreetMap data are in geographic coordinates the length is in degree, but I prefer to have the length of a way in meter. Thankfully PostGIS provides ST_Distance_Spheroid to calculate lengths on the (earth) spheorid:

UPDATE ways
SET length = ST_Length_Spheroid(the_geom,'SPHEROID["WGS 84",6378137,298.257223563]')::DOUBLE PRECISION;

Next step is to filter the one-ways paths and set the "reverse_cost" to -1 to force pgRouting to choose these ways only in forward direction. Unfortunately only the name and the osm_id are imported from the OpenStreetMap raw data to the database. Thus I needed a list of osm_ids that represents one-ways. Since I also used a PostgreSQL database to store the GRASS attribute tables while working with OpenStreetMap data, I could easily get a comma-separated list of osm_ids with the following query:

psql -A -t -R ',' -d grass -U grass -c "SELECT osm_id FROM osm_net WHERE oneway = 'yes';" > oneways.sql

Then I updated the "reverse_cost" column:

UPDATE ways SET reverse_cost = -1
WHERE osm_id IN(55909972,55909973,28977387, ... ,102306365,27158132)

Finally I could start with undirected and directed queries and compare the results.

-- Undirected shortest path query with astar
SELECT * FROM astar_sp( 'ways',3365,3412)
-- Directed shortest path query with astar
SELECT * FROM astar_sp_directed( 'ways',3365,3412,true,true)

Verify both calculated shortest paths from node 3365 to node 3412 in QGIS:

The red path is the undirected query and ignores apparently one-ways, contrary to the directed green path that follows one-way directions.

Until now I was considering only length and direction of a way. To get a more life-like model I wanted to take the road classes into account, too. I updated table "public.classes" and estimated for each road class an average speed in m/s in column "cost":

psql -A -d osm_routing -U openstreet -c "SELECT id,type_id,TRIM(trailing FROM name) AS name,cost FROM classes LIMIT 8;"
id|type_id|name|cost
101|1|motorway|33.3333
102|1|motorway_link|27.7778
103|1|motorway_junction|27.7778
104|1|trunk|27.7778
105|1|trunk_link|27.7778
106|1|primary|13.8889
107|1|primary_link|13.8889
108|1|secondary|11.1111
(8 rows)

With these assumed costs for each road class I updated the costs in table "public.ways" and divided the length by the average speed and got the time in seconds that is necessary to pass a certain way.

-- Update to_cost column
UPDATE ways
SET to_cost=(length/(SELECT "cost" FROM classes WHERE id=class_id))::DOUBLE PRECISION
-- Update reverse_cost columns
UPDATE ways
SET reverse_cost=(length/(SELECT "cost" FROM classes WHERE id=class_id))::DOUBLE PRECISION
WHERE reverse_cost <> -1;

Now everything is set up and it needs nothing but another rainy weekend to compare pgRouting results with GRASS GIS results.

Last but not least I'd like to share yet another way how to display results of PostGIS queries without using an intermediate (Shape-)file, namely using the very versatile Virtual Format from the OGR library. Create a text file with the following lines and open it as a vector layer in QGIS or any other OGR backed application:

<OGRVRTDataSource>
  <OGRVRTLayer name="shortest_path">
    <SrcDataSource>
      PG:dbname='osm_routing' user='openstreet'
    </SrcDataSource> 
    <SrcSQL>
      SELECT * FROM astar_sp_directed( 'ways',3365,3414,true,true)
    </SrcSQL>
  </OGRVRTLayer>
</OGRVRTDataSource>

GPX download available on openstreetmap.la

I'm pleased to announce minor updates on openstreetmap.la: OpenStreetMap points of interest in Laos are now available for download in the popular GPS Exchange Format (GPX) and the roads Shapefile has been improved.

I considered providing GPX files as an additional download format since the launch of openstreetmap.la, but I didn't implemented it till last week. Contrary to the KMZ format, that always requires (e.g. PHP) scripting if you want customized styles, my goal was to create the GPX file without any scripting.

As often there are several ways how to solve a problem, I'd like to share how I did it, assuming the OpenStreetMap are already in a PostGIS database. As a (probably) daily user of OGR, I knew that OGR supports PostGIS as well as GPX, so it was the obvious way to have a deeper look at OGR's GPX driver.

Since GPX uses a fixed schema the attributes from the database has to be renamed according to the GPX schema. It was obvious how to name the attributes to get the <name>, <cmt> etc. elements. But from the OGR documentation I didn't understand how I have to name the link attribute, since the <link> element is nested and the GPX schema allows as many <link> elements as you want. I took a closer look at the code to see that OGR requires the first link named "link1_href", the second "link2_href" etc. Furthermore I wanted different symbols for the different OpenStreetMap features like restaurant, hotel etc. that required a nested CASE WHEN ... ELSE statement. Since the SQL statement was getting quite impressive, I decided to create a new database view called gpx_poi, instead of using the -sql OGR option or creating a new virtual OGR format.

CREATE OR REPLACE VIEW gpx_poi AS
SELECT "way",
CASE WHEN "name" IS NOT NULL THEN "name"
  ELSE CASE WHEN "amenity" IS NOT NULL THEN INITCAP(REPLACE("amenity",'_',' '))
    ELSE CASE WHEN "tourism" IS NOT NULL THEN INITCAP(REPLACE("tourism",'_',' '))
      ELSE ''
    END
  END
END AS "name",
CASE WHEN "amenity" IS NOT NULL THEN INITCAP(REPLACE("amenity",'_',' '))
  ELSE CASE WHEN "tourism" IS NOT NULL THEN INITCAP(REPLACE("tourism",'_',' '))
    ELSE ''
  END
END AS "cmt",
CASE WHEN "amenity" IS NOT NULL THEN INITCAP(REPLACE("amenity",'_',' '))
  ELSE CASE WHEN "tourism" IS NOT NULL THEN INITCAP(REPLACE("tourism",'_',' '))
    ELSE ''
  END
END AS "desc",
"link1_href",
"link1_text",
'text/html'::TEXT AS "link1_type",
CASE WHEN "tourism" = 'hotel' THEN 'Lodging'
  ELSE CASE WHEN "amenity" = 'post_office' THEN 'Post Office'
    ELSE CASE WHEN "amenity" = 'telephone' THEN 'Telephone'
      ELSE CASE WHEN "amenity" = 'toilets' THEN 'Restrooms'
        ELSE CASE WHEN "amenity" = 'school' THEN 'School'
-- etc. for all other symbols
        END
      END
    END
  END
END AS "sym"
FROM (SELECT "way","osm_id","name","amenity",
"tourism","natural","man_made","place",
'http://www.openstreetmap.org/browse/node/' || "osm_id" AS "link1_href",
'Node: ' || "osm_id" AS "link1_text"
FROM public.osm_en_point
UNION
SELECT ST_Centroid(way) AS "way","osm_id","name","amenity",
"tourism","natural","man_made","place",
'http://www.openstreetmap.org/browse/way/' || "osm_id" AS "link1_href",
'Way: ' || "osm_id" AS "link1_text"
FROM public.osm_en_polygon)
AS union_table

Having this gpx_poi database view, it's an easy ogr2ogr conversion:

ogr2ogr -f GPX -s_srs EPSG:900913 -t_srs EPSG:4326 pointsofinterest.gpx PG:"dbname='openstreet' user='openstreet'" gpx_poi

Just to make sure the GPX output is correct, I did also an XML validation of the output file:

xmlstarlet val --xsd http://www.topografix.com/GPX/1/1/gpx.xsd --err pointsofinterest.gpx

And finally the result, a screenshot from the waypoints on a Garmin device:

The second update on openstreetmap.la concerns the roads Shapefile. Last week when I worked on the routing in GRASS GIS I realized that the road reference number is missing in the roads Shapefile. Meanwhile the "ref" attribute is included.