Saturday, July 30, 2011

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

Thursday, July 28, 2011

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.


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


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

Friday, July 15, 2011

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:
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;"
(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
SET to_cost=(length/(SELECT "cost" FROM classes WHERE id=class_id))::DOUBLE PRECISION
-- Update reverse_cost columns
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:
  <OGRVRTLayer name="shortest_path">
      PG:dbname='osm_routing' user='openstreet'
      SELECT * FROM astar_sp_directed( 'ways',3365,3414,true,true)

Monday, July 4, 2011

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 ""
cd shapelib-1.3.0b2/contrib/
wget ""
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):
< 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
< 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 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):
<     identifyKey = DBFGetFieldIndex( iDBF, "OSM_ID" );
<     if( identifyKey >= 0 )
<         goto gotkey;
<     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  
TOURISM: guest_house
NAME: Annivong 2
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
NAME: restaurant
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:
Screenshot QGIS Attribute 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.

Saturday, July 2, 2011

GPX download available on

I'm pleased to announce minor updates on 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, 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.
SELECT "way",
  ELSE CASE WHEN "amenity" IS NOT NULL THEN INITCAP(REPLACE("amenity",'_',' '))
    ELSE CASE WHEN "tourism" IS NOT NULL THEN INITCAP(REPLACE("tourism",'_',' '))
      ELSE ''
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 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 AS "desc",
'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 AS "sym"
FROM (SELECT "way","osm_id","name","amenity",
'' || "osm_id" AS "link1_href",
'Node: ' || "osm_id" AS "link1_text"
FROM public.osm_en_point
SELECT ST_Centroid(way) AS "way","osm_id","name","amenity",
'' || "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 --err pointsofinterest.gpx
And finally the result, a screenshot from the waypoints on a Garmin device:

The second update on 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.