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

Search bar on openstreetmap.la

I'm pleased to announce a new feature on openstreetmap.la: a search bar to find places and points of interests like hotels, hospitals, restaurants or convenience shops.

To facilitate the search, the search bar implements auto-completion: Start typing the name of any place or point of interest, and a drop-down menu suggests available places or points containing this name. A small icon in front of the name indicates if it's a village, restaurant, hotel etc. The icons corresponds in most (!) cases to the map symbols. Villages, towns and cities are the exception, since places have a label but no icon on the map. The icons for places used in the drop-down menu look as follows:

Village

Town

City

See the following screenshots:

Search bar with autocompletion

Works also in Lao

A map marker highlights the found place:

Map marker

As always: feedback is highly appreciated!

Open Source Routing Machine: Yet another routing engine

Initially I only wanted to some routing on OpenStreetMap data with GRASS GIS. Then I also tried pgRouting and later SpatiaLite. Using pgRouting as backend, I set up a simple web application following the well explained pgRouting workshop. But I wasn't really satisfied since the routing starts only at the closest node, that's basically less a problem in the city than on the countryside with less junctions (and thus less nodes).

Finally I turned away from GIS to specialized routing software for OpenStreetMap, namely to Routino and to the Open Source Routing Machine, short OSRM. I was that impressed by the speed of OSRM, that I decided to have a deeper look at this engine:

Pros

• Speed! OSRM is about speed and it claims to be capable to handle continental sized networks using memory efficient third party libraries like google-sparsehash. It was a matter of seconds to preprocess the 6.2 MB Laos PBF file.
• Web-oriented: includes a server that runs on a configurable port and thus it is very easy to integrate in any web application.
• Since it's developed for OpenStreetMap, it considers out of the box OpenStreetMap tags like oneways etc.

Cons

• Configuration is not trivial, since speed and highway categories are hard-coded.
• Installation with Scons was very laborious on the server, since I had some dependencies in my home directory and it was necessary to fiddle the SConstruct file.
• Missing GeoJSON support. But since it needs anyway a proxy script to use OSRM in Ajax requests, it was a quick hack to reformat the homebrew JSON output to proper GeoJSON.

I adapted the simple web frontend I developed for pgRouting to work with OSRM as backend and put everything online. I'm pleased to present the current proof-of-concept prototype:

http://www.openstreetmap.la/dev/routing/

As soon as time allows I'll integrate the routing on the openstreetmap.la main page.

Power lines from Landsat imagery

Hydropower is the big deal in Laos at the moment and therefore power lines are required. Some weeks ago on an offroad ride we crossed the relatively new Nam Ngum 2 power transmission line. Realizing the long straight and wide forest aisle, I knew that it must be possible to identify the line on Landsat imagery to map it on OpenStreetMap. Landsat imagery is in the public domain and thus suitable for mapping.

I downloaded the Landsat imagery from January and February 2011 for this region already earlier. Therefore I could start mapping the power line immediately, starting from the point where the power line crossed our way and following the forest aisle on Landsat. It's probably the first mapped power line on OpenStreetMap in Laos.

Nam Ngum 2 dam and power lines on Landsat

There is the Nam Ngum 2 reservoir and its dam visible in the middle at the top border. In the lower left corner is the older Nam Ngum 1 reservoir. The power line going south is clearly identifiable in the middle section. Furthermore two access roads to the dam are identifiable, one on the north shore of Nam Ngum 1 and one leading east from the dam. Just downstream of the new dam there is the confluence of Nam Bak.

Just for the visual appearance I tweaked the image to get rid of the stripes using GDAL, GIMP and OSSIM. OSSIM is as fas as I know still the only FOSS that implements histogram stretch based on the standard deviation.

The power lines in natura

Improved Polygon Capturing 1.0 is released

I'm pleased to announce version 1.0 of the Improved Polygon Capturing plugin for QGIS!

Until now I didn't want to give version 1.0 to this plugin since one important feature was still missing: After digitizing a new geometry the feature form was not opened like it is the usual behavior of the standard editing tools. But today I discovered very coincidentally the openFeatureForm method in the QgisInterface class. It wasn't a big deal to integrate it and now the feature form is opened and the plugin gets to version 1.0.
Feedback is highly appreciated!

From the help file (last update 8.8.2011):

Improved Polygon Capturing is a QGIS Python plugin, that allows to digitize new polygons or lines with predefined edge lengths. Point layers are not handled.

The plugin is available in the QGIS contributed repository. Please see the official QGIS manuals to get further information about plugin repositories.

How to use

The plugin adds a new icon and a spin box to the digitizing toolbar. The icons are derived from the gis icon theme and look as follows:

Editing a polygon layer:	Editing a line layer:

The digitizing toolbar with added plugin:

Select a polygon or a line layer and toggle editing. After setting the first vertex the next vertices are added in the direction of the mouse pointer with the predefined length set in the spin box. Similar to the standard editing tools left mouse clicks add a new vertex while right mouse clicks finish the geometry.
After finishing a new geometry the feature form opens and attributes can be entered.
While digitizing new vertices the snapping properties from the current project settings are considered as well as the avoiding intersection properties.
If the distance in the spin box is set to 0 (zero), new vertices are set at the mouse position.

Caveats

The plugin calculates the distance in plain trigonometry. Thus it is not recommended to use it in unprojected systems like EPSG:4326.

History

Improved polygon capturing plugin has been developed in June 2010 for a land management and registration project in the Lao PDR. It is used to digitize land parcels with known edge lengths from high-resolution satellite images. The parcel edges have been measured a priori in the field.

• Version 0.8 (June 2010): first published version
• Version 0.9 (February 2011): added support for polyline layers
  and bug fixing
• Version 1.0 (July 2011): feature form opens in editing
  mode after finishing a new feature
• Version 1.0.1 (August 2011): plugin gets an icon in the plugin manager and rubber band line width and color settings are considered

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>

Comparison of Shapefiles with shpdiff

For my important coding projects I'm using a local Subversion repository and a file comparison utility like diff (provided by the IDE) to keep track of my changes.

Some time ago I started to think about a similar versioning system for Shapefiles, since in one of our projects we have GIS staff concurrently tracing and updating base data from satellite imagery. Finally I found the shpdiff utility, which is exactly what I looked for. These were my build steps following a thread from the ShapeLib mailing list.

Download and extract the latest shapelib version and the shpdiff.c file. Probably there are also other places where you can find the shpdiff.c file.

wget "http://download.osgeo.org/shapelib/shapelib-1.3.0b2.zip"
unzip shapelib-1.3.0b2.zip
cd shapelib-1.3.0b2/contrib/
wget "ftp://ftp.soc.soton.ac.uk/pub/pwc101/slackware/slackbuilds/academic/shp2text/shpdiff.c"

Next step is to build the main utilities and the library

cd shapelib-1.3.0b2/
make all
make lib

Then I wanted to build the community contributed tools. Before compiling I made the following changes to the Makefile to include the shpdiff utility (diff output):

18c18
< all: shpdxf shpproj dbfinfo shpcentrd shpdata shpwkb dbfinfo dbfcat shpinfo shpfix shpcat Shape_PointInPoly shpsort shpdiff
---
> all: shpdxf shpproj dbfinfo shpcentrd shpdata shpwkb dbfinfo dbfcat shpinfo shpfix shpcat Shape_PointInPoly shpsort
37,39d36
< 
< shpdiff: shpdiff.c $(SHPOBJ)
<  $(CC) $(CFLAGS) shpdiff.c ${SHPOBJ} $(LINKOPT) $(GEOOBJ) -o shpdiff

finally I built also these tools with

make all

It seems that the changes in shputil.c and shpgeo.c described in the mentioned thread are no longer necessary. To test shpdiff I used the amenities Shapefile from openstreetmap.la. I started without any changes:

./shpdiff amenities.shp amenities_modified.shp
Original Shapefile Type: Point, 1392 shapes 1392 database records
Comparison Shapefile Type: Point, 1392 shapes 1392 database records

NOTE: Using column NAME to identify records

It's an important note by shpdiff that it compares features based on their "NAME" attribute (and then "STREET", "TOWN" etc.). Since OpenStreetMap features have an unique identifier, I wanted to compare the files based on the "OSM_ID" attribute. These are my changes in shpdiff.c on lines 173 and 208 (again diff output):

173,175d172
<     identifyKey = DBFGetFieldIndex( iDBF, "OSM_ID" );
<     if( identifyKey >= 0 )
<         goto gotkey;
208c205
<     if(identifyKey >= 0)
---
>     if(identifyKey)

If you have another identifier attribute you want to use for comparison, just change the name and rebuild it:

make shpdiff

Then I started to edit one of the files. First I deleted a feature, shpdiff outputs the following:

Original Shapefile Type: Point, 1392 shapes 1392 database records
Comparison Shapefile Type: Point, 1391 shapes 1391 database records

NOTE: Using column OSM_ID to identify records

Record 948: deleted from original
--------------------------------
OSM_ID: 513607024  
AMENITY:                                     
TOURISM: guest_house
NAME: Annivong 2
NAME_LO:                                                             
NAME_EN: Annivong 2

The output after adding a new feature (new OSM ID should be negative, right?):

Original Shapefile Type: Point, 1392 shapes 1392 database records
Comparison Shapefile Type: Point, 1393 shapes 1393 database records

NOTE: Using column OSM_ID to identify records

New record 1392 found
-------------------------------
OSM_ID: -1
AMENITY: restaurant
TOURISM:
NAME: restaurant
NAME_LO:
NAME_EN:

Move a feature:

Original Shapefile Type: Point, 1392 shapes 1392 database records
Comparison Shapefile Type: Point, 1392 shapes 1392 database records

NOTE: Using column OSM_ID to identify records

Record 832:shape change

Change an existing field value for feature 387, delete a value for feature 667 and add a value for feature 717

Original Shapefile Type: Point, 1392 shapes 1392 database records
Comparison Shapefile Type: Point, 1392 shapes 1392 database records

NOTE: Using column OSM_ID to identify records

Record 387:
NAME_EN: Mittapharb Lao Barbecue >>> Mittaphab Lao Barbecue

Record 667:
NAME_EN: Ho Phra Keo >>> 

Record 717:
NAME_EN:  >>> Wat Kaognot

The output is quite self-explanatory and you can find the records in the attribute table with the corresponding number as long as you don't sort the table:

For the sake of completeness, I should also mention the PostGIS Versioning scripts and QGIS plugin from Kappasys. It's not only a comparison tool but a whole versioning system and looks very promising and sophisticated. But since Shapefile is still the dominant GIS format in my working environment I haven't yet tried it.