OrbView-3 imagery freely available from USGS

geo-spatial.org's mailing list had a link to this little gem: High Resolution OrbView - 3 Images Now Available from USGS.

The resolution itself is not spectacular when compared to the Bing imagery displayed in OSM's Potlatch and the coverage of Romania's territory is rather sparse, but the license lets you do whatever you want with the imagery (public domain).

A preview of the OSM street network on top of OrbView-3 imagery for Bucharest, Aviatiei:

geo-spatial.org are working on a project to make the processed data available through their WMS service.

Line of sight analysis using Grass GIS and SRTM elevation data

Grass GIS has built-in support for line of sight analysis (the command is called r.los).

I'm revisiting an older post (Mountain Hiking Tracks in Făgăraș) that I've ended with a line of sight diagram I've built by hand based on the description given by a fellow hiker who described to us what we were seeing (peak names). This time I've used Grass GIS to build the line of sight raster from Vânătoarea lui Buteanu peak with a 20 Km radius. The background is based on the freely available Landsat 7 image data.

Note there are NULL values (failures to read altitudes in certain raster cells) in the SRTM data (represented in black in the picture below) so the line of sight analysis above was carried on an interpolated version of the data.

As far as I can tell, the diagram is fairly accurate. I'm not sure about "Turnul Paltinului" and the small ridge following it, but the rest of the features fall into place pretty nicely. This is a photo taken from Vânătoarea lui Buteanu towards Negoiu Peak:

And this is a mash-up of the line of sight diagram over-imposed on the photograph:

Using OpenStreetMap osm2pgsql data in Grass GIS

These are my notes on using the OSM data as imported by osm2pgsql in a PostgreSQL database with Grass GIS.

While Grass is a very powerful tool, it's Unixy feel is a little disconcerting. If you're a programmer like me, it helps to realize the Grass command line shell is in fact a preconfigured bash session:

GRASS 6.4.1 (OSM_900913):~ > ps ax | grep "$$" | grep -v grep
98804 s000  S      0:00.06 /bin/bash

The things you get from Grass are a set of binaries that interact with the data you have (organized in locations that are part of a mapset stored inside your filesystem). The other thing that is somewhat unexpected for a PostGIS user is that, within a Grass location, all the data has to share a single projection.

Creating a location for your OSM data

The first step required when using OSM data is to create a Grass location for the OSM data. I've used the same projection osm2pgsql uses when importing, the "fake" projection with SRID 900913. Grass doesn't know it, so you can define it as a custom PROJ4 string based on the definition you get from PostGIS:

osm=# select proj4text from spatial_ref_sys where srid = 900913;
                                                    proj4text                                                     
------------------------------------------------------------------------------------------------------------------
 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs
(1 row)

Here are the relevant screen shots when defining the location using Grass Python WXGUI:

Creating temporary tables in OSM

To reduce the size of the data set, I've created a couple of tables in PostGIS by selecting only those features I planned to use for my map by spatially selecting those features within a certain administrative boundary:

create table bucharest_s1_streets as
select * from planet_osm_line where st_contains((select way from planet_osm_polygon where boundary = 'administrative' and name ='Sector 1' limit 1), way)
and highway is not null;

create table bucharest_s1_pois as
select * from planet_osm_point where st_contains((select way from planet_osm_polygon where boundary = 'administrative' and name ='Sector 1' limit 1), way);

I then defined these tables in PostGIS' geometry_tables (otherwise Grass will refuse to import them):

insert into geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type) values('','public','bucharest_s1_pois','way',2,900913,'POINT');
insert into geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type) values('','public','bucharest_s1_streets','way',2,900913,'LINESTRING');

Importing the data in Grass

For some reason, my Grass setup fails to set the environment variables below, so I have to re-export them:

GRASS 6.4.1 (OSM_900913):~ > export GISDBASE=/Users/diciu/grassdata
GRASS 6.4.1 (OSM_900913):~ > export LOCATION_NAME=OSM_900913
GRASS 6.4.1 (OSM_900913):~ > export MAPSET=PERMANENT

The PostGIS database parameters have to be defined next:

export HOST=localhost
export DB=osm
export USER=gisuser
export PASS=somepassword
export DSN="PG:host=${HOST} dbname=${DB} user=${USER} password=${PASS}"

You can now run the v.ogr.in commands to import your data (I'm importing roads and pharmacies):

v.in.ogr dsn="${DSN}" layer=bucharest_s1_streets out=streets
v.in.ogr dsn="${DSN}" layer=bucharest_s1_pois out=pharmacy where="amenity='pharmacy'"

If everything went well, not you have two vector layers in your location:

GRASS 6.4.1 (OSM_900913):~ > g.list type=vect
----------------------------------------------
vector files available in mapset :
pharmacy                  streets
----------------------------------------------

Computing distance to pharmacies

The Vector networking tutorial describes a number of very nice examples of processing tools available in Grass. I am using v.net.iso below to compute the distance from pharmacies:

v.distance -p from=pharmacy to=streets output=pharmacy_conn_streets upload=dist column=dist
v.patch in=pharmacy,streets,pharmacy_conn_streets out=pharmacy_net
v.clean in=pharmacy_net out=pharmacy_net_clean tool=snap,break thresh=1
v.net.iso input=pharmacy_net_clean output=pharmacy_iso ccats=1-100 costs=200,500,1000 nlayer=1

The command's output are a number of vector layers, pharmacy_iso being the one that stores our road network categorized by distance to pharmacies (less than 200 meters, between 200 and 500 meters and so on).

Rendering a map using the data

It's much easier and the results are probably better if you use QGIS instead of Grass' UI. There's a plugin for QGIS that asks you to specify the Grass mapset, location and loads a given layer. As an added bonus, QGIS is able to re-project on the fly so you can mix layers in different projections.

Here's a map of pharmacies in my neighborhood as rendered by QGIS' export to PNG feature:

Using pgRouting on an OSM database imported by osm2pgsql

If you've tried to use pgRouting on an OSM database imported by osm2pgsql you've probably noticed it doesn't work.

It doesn't work because one of the requirements for pgRouting is that streets are split at intersections, meaning it expects each osm_id from planet_osm_line to only contain street segments.
However, in planet_osm_line, a street is represented by a LINESTRING, most of the times the intersection with other streets being just another POINT in the street's LINESTRING.

A visual representation of the long explanation above:







The good news is that it's fairly easy to split the OSM streets into constituent segments using PLPGSQL.
One way to achieve that is the one listed below (note the filtering on "Sector 1" - I am restricting table entries within the administrative boundary of the first sector of Bucharest because things get really slow on the entire OSM road network):

drop table if exists network;
create table network(gid serial, osm_id integer, name varchar, the_geom geometry, source integer, target integer, length float);


CREATE OR REPLACE FUNCTION compute_network() RETURNS text as $$
DECLARE
streetRecord record;
wayRecord record;
pointCount integer;
pointIndex integer;
geomFragment record;
BEGIN
-- for each street
FOR streetRecord in select way, osm_id, name from planet_osm_line where highway is not null and st_contains((select way from planet_osm_polygon where boundary = 'administrative' and name like 'Sector 1' limit 1), way) LOOP
 -- for each street in the region of interest
 SELECT * from planet_osm_ways where id = streetRecord.osm_id into wayRecord; 
 FOR pointIndex in array_lower(wayRecord.nodes, 1)..array_upper(wayRecord.nodes,1)-1 LOOP
  RAISE NOTICE 'Inserting name % source %, target %', streetRecord.name, wayRecord.nodes[pointIndex], wayRecord.nodes[pointIndex+1];
  select st_makeline(st_pointn(streetRecord.way, pointIndex), st_pointn(streetRecord.way, pointIndex+1)) as way into geomFragment;
  insert into network(osm_id, name, the_geom, source, target, length) values(streetRecord.osm_id, streetRecord.name, geomFragment.way, wayRecord.nodes[pointIndex], wayRecord.nodes[pointIndex+1], st_length(geomFragment.way));
 END LOOP;
END LOOP;
return 'Done';
END;
$$ LANGUAGE 'plpgsql';

select * from compute_network();
select assign_vertex_id('network', 1, 'the_geom', 'gid');

The PLPGSQL procedure works by creating a street segment for each pair of nodes from planet_osm_ways; the geometry is created by simply creating a new line from the points corresponding to the nodes (see the st_makeline call).

Once assign_vertex_id has been run, you can start using shortest_path:

osm=# select n.name, p.edge_id, o.highway from network n, (select edge_id from shortest_path('
osm'# SELECT gid as id,
osm'# source,
osm'# target,
osm'# length as cost
osm'# FROM network',
osm(# (select source from network where name='Intrarea Jugului' limit 1),(select source from network where name='Strada Sofia' limit 1),false,false)) p,
osm-# planet_osm_line o
osm-# where p.edge_id = n.gid and n.osm_id = o.osm_id;
             name             | edge_id |   highway   
------------------------------+---------+-------------
 Intrarea Jugului             |    2330 | residential
 Strada Elena Caragiani       |     808 | residential
 Strada Borșa                 |    2016 | residential
 Strada Borșa                 |    2015 | residential
 Strada Borșa                 |    2014 | residential
 Strada Borșa                 |    2013 | residential
                              |    1992 | residential
 Șoseaua Pipera               |    8098 | secondary
 Strada Nicolae Caramfil      |    1828 | secondary
 Strada Nicolae Caramfil      |    1827 | secondary
 Strada Nicolae Caramfil      |    1826 | secondary
 Strada Nicolae Caramfil      |    1825 | secondary
 Strada Nicolae Caramfil      |    1824 | secondary
 Strada Nicolae Caramfil      |    1823 | secondary
 Strada Nicolae Caramfil      |    1822 | secondary
 Strada Nicolae Caramfil      |    1821 | secondary
 Bulevardul Beijing           |    1035 | secondary
 Bulevardul Beijing           |    1036 | secondary
 Bulevardul Beijing           |    7883 | secondary
 Bulevardul Beijing           |    7766 | secondary
 Bulevardul Beijing           |    1474 | residential
 Bulevardul Beijing           |    1473 | residential
 Bulevardul Beijing           |    1472 | residential
 Bulevardul Beijing           |    1471 | residential
 Bulevardul Beijing           |    1470 | residential
 Bulevardul Beijing           |    1469 | residential
 Bulevardul Beijing           |    1468 | residential
 Bulevardul Beijing           |    1467 | residential
 Piața Charles de Gaulle      |     859 | residential
 Piața Charles de Gaulle      |     858 | residential
 Piața Charles de Gaulle      |    8918 | residential
 Piața Charles de Gaulle      |    8917 | residential
 Piața Charles de Gaulle      |    8916 | residential
 Piața Charles de Gaulle      |    8915 | residential
 Piața Charles de Gaulle      |    8914 | residential
 Piața Charles de Gaulle      |    8913 | residential
 Piața Charles de Gaulle      |    8912 | residential
 Calea Dorobanților           |     264 | primary
                              |    9999 | footway
 Calea Dorobanților           |    8116 | primary
 Calea Dorobanților           |    8117 | primary
 Calea Dorobanților           |    8118 | primary
 Calea Dorobanților           |    8119 | primary
 Strada Emil Pangrati         |    3489 | residential
 Strada Emil Pangrati         |    3412 | residential
 Strada Emil Pangrati         |    3413 | residential
 Strada Emil Pangrati         |    3414 | residential
 Cpt. Av. Demetriade Gheorghe |    3446 | residential
 Strada Andrei Muresanu       |    8086 | residential
(49 rows)

Time: 77.807 ms

If you render the above (I used QGIS), the (start of) the route looks like this:



Note how pgRouting has found multiple edges (2016, 2015, 2014, 2013) alongside "Strada Borșa"; that's because we split the street using the PLPGSQL procedure described at the beginning of the post:

osm=# select osm_id, name from network where gid in (2016, 2015, 2014, 2013);
  osm_id  |     name     
----------+--------------
 23730998 | Strada Borșa
 23730998 | Strada Borșa
 23730998 | Strada Borșa
 23730998 | Strada Borșa
(4 rows)

Installing the Fast SQL Layer plugin on Mac OS X 10.7

If you're running QGis version 1.7 (trunk) from the KyngChaos site and you want to install the Fast SQL Layer plugin, you need to solve some dependencies first.

The paths below may not be the same on your machine if you've not installed GEOS and PROJ4 from the KyngChaos packages.

sudo easy_install-2.6 Pygments
sudo C_INCLUDE_PATH=/usr/local/pgsql/include/ easy_install-2.6  psycopg2
sudo LIBRARY_PATH=/Library/Frameworks/GEOS.framework/Versions/3/unix/lib/:/Library/Frameworks/PROJ.framework/Versions/4/unix/lib/ C_INCLUDE_PATH=/Library/Frameworks/GEOS.framework/Versions/3/unix/include/:/Library/Frameworks/PROJ.framework/Versions/4/Headers/ easy_install-2.6  pyspatialite

If the installation has been successful, you should not get an error when importing the modules, i.e.:

cristi:~ diciu$ python2.6
Python 2.6.6 (r266:84292, Jun 16 2011, 16:59:16) 
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pygments, psycopg2, pyspatialite
>>> 

Once you've installed the dependencies, you can go on installing the plugin using the QGIS GUI.

Using pgRouting to compute routes on OSM hiking trails

I've experimented with using pgRouting's shortest_path function to compute routes using the hiking trails defined in the OpenStreetMap data set.

Building the topology

In the OSM model as created by osm2pgsql, the relations table models a hiking trail based on a set of members (each member a LINESTRING geometry from the planet_osm_line table).

The problem is you cannot build a topology based solely on the start and end points of each trail segment, because more often then not, trail segment intersect mid-way:




So I used the relations table to build the topology where each trail segment is split based on the POINTS that are part of its LINESTRING geometry. This means that every point from the trail segment becomes a node in our topology (this is very inefficient, but I could not find a faster way of not missing mid-trail intersections).

I've built the topology using the following SQL code:

drop table if exists trail_network;

create table trail_network(gid serial, osm_id integer, name varchar, symbol varchar, the_geom geometry, source integer, target integer, length float);



CREATE OR REPLACE FUNCTION compute_trail_network() RETURNS text as $$

DECLARE

relationRecord record;

partIndex integer;

lineRecord record;

wayRecord record;

geomFragment record;

BEGIN

	FOR relationRecord in select id, parts, tags from planet_osm_rels where tags @> '{osmc:symbol}' and tags @> '{hiking}' LOOP

		FOR partIndex in array_lower(relationRecord.parts, 1)..array_upper(relationRecord.parts, 1) LOOP

			SELECT name, highway, way from planet_osm_line where osm_id = relationRecord.parts[partIndex] into lineRecord;

			SELECT * from planet_osm_ways where id = relationRecord.parts[partIndex] into wayRecord;

			FOR pointIndex in 1..ST_NPoints(lineRecord.way)-1 LOOP

				select st_makeline(st_pointn(lineRecord.way, pointIndex), st_pointn(lineRecord.way, pointIndex+1)) as way into geomFragment;

				INSERT into trail_network(osm_id, name, symbol, the_geom, length) values(relationRecord.parts[partIndex], substring(array_to_string(relationRecord.tags, '#'), '#name#([^#]*)#?$'), substring(array_to_string(relationRecord.tags, '#'), '(red_cross|blue_cross|red_triangle|blue_triangle|red_dot|blue_dot|red_stripe|blue_stripe|yellow_stripe|yellow_triangle|yellow_cross)'), geomFragment.way, st_length(geomFragment.way));

			END LOOP;

		END LOOP;

	END LOOP;

return 'Done';

END;

$$ LANGUAGE 'plpgsql';



select * from compute_trail_network();

Once the topology has been created, I called assign_vertex_id() to compute source and target values (for some reason, probe_geometry_columns does not work on my machine so I had to fill in my topology table's geometry column data by hand):

insert into geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type) values('', 'public', 'trail_network', 'the_geom', 2, 900913, 'LINESTRING');



SELECT assign_vertex_id('trail_network', 2, 'the_geom', 'gid');

Computing a route

The first thing required is finding topology nodes close to the POI's that represent departure and arrival points.
In this example, we will compute a route between Cabana Piscul Negru and Negoiu:

osm=# select t.gid, t.source, t.target, t.osm_id, t.name, st_distance(p.way, t.the_geom) from trail_network t, planet_osm_point p where st_buffer(p.way, 100) && t.the_geom and p.name='Cabana Piscul Negru' order by st_distance(p.way, t.the_geom) asc limit 5;

  gid  | source | target |  osm_id  |              name              |   st_distance    

-------+--------+--------+----------+--------------------------------+------------------

 41354 |  17684 |  17685 | 73964460 | Piscu Negru - Refugiul Călțun  | 94.9336094118476

 41532 |  17684 |  17685 | 73964460 | Piscul Negru - Șaua Paltinului | 94.9336094118476

(2 rows)

there are two edges in close proximity to Piscul Negru (94 meters away). They share source and target nodes because they are in fact over a segment shared between the two different trails.

We repeat the query to find out the node for Negoiu:

osm=# select t.gid, t.source, t.target, t.osm_id, t.name, st_distance(p.way, t.the_geom) from trail_network t, planet_osm_point p where st_buffer(p.way, 100) && t.the_geom and p.name='Negoiu' order by st_distance(p.way, t.the_geom) asc limit 5;

  gid  | source | target |  osm_id  |      name      |   st_distance    

-------+--------+--------+----------+----------------+------------------

 42688 |  18832 |  18833 | 27182093 | Traseu creasta | 2.03276684697733

 42687 |  18831 |  18832 | 27182093 | Traseu creasta | 2.18879960182101

 42686 |  18830 |  18831 | 27182093 | Traseu creasta | 13.9052808592012

 42685 |  18829 |  18830 | 27182093 | Traseu creasta | 35.2142772300472

 42684 |  18828 |  18829 | 27182093 | Traseu creasta |  63.933925517158

(5 rows)

Now that we've got the start (17684) and end (18832) nodes, we can ask pgRouting to compute the route:

osm=# select tn.name, p.edge_id, tn.symbol from trail_network tn, (select edge_id from shortest_path('

osm'# SELECT gid as id,

osm'# source,

osm'# target,

osm'# length as cost

osm'# FROM trail_network',

osm(# 17684, 18832,false,false)) p

osm-# where p.edge_id = tn.gid;

                name                 | edge_id |    symbol     

-------------------------------------+---------+---------------

 Piscu Negru - Refugiul Călțun       |   41343 | blue_triangle

 Piscu Negru - Refugiul Călțun       |   41341 | blue_triangle

 Piscu Negru - Refugiul Călțun       |   41338 | blue_triangle

 Piscu Negru - Refugiul Călțun       |   41335 | blue_triangle

[..]

If we retrieve a summary of the involved trails, we get:

osm=# select max(tn.name), max(tn.symbol) from trail_network tn, (select edge_id from shortest_path('

osm'# SELECT gid as id,

osm'# source,

osm'# target,

osm'# length as cost

osm'# FROM trail_network',

osm(# 17684, 18832,false,false)) p

osm-# where p.edge_id = tn.gid group by tn.name;

                 max                 |      max      

-------------------------------------+---------------

 Piscu Negru - Refugiul Călțun       | blue_triangle

 Lacul Călțun - Tunel Transfăgărășan | blue_cross

 Traseu creasta                      | red_stripe

(3 rows)

While the choice of the blue cross trail might look weird, it happens because in this particular section, the blue cross trail shares the same path as the blue triangle trail. From the routing engine's point of view, they are perfectly equivalent but a human might argue that the correct trail sequencing would be "Piscu Negru - Refugiul Caltun" following the blue triangle and then the ridge trail following the red stripe up to Negoiu.

You can also use the driving distance functions from pgRouting.
This is a colormap showing a classification of reachable places starting from Cabana Piscul Negru:

Piatra Craiului mashup

A mashup of the OpenStreetMap hiking trails displayed in Google Earth:



I've initially experimented with displaying the actual OSM tiles over Google Earth (that automatically maps them using the terrain's altitude profile) but I obtained mixed results - for one, the tiles overlapped the satellite imagery displayed by GE and, most annoyingly, I didn't manage to find a way to correctly match GE's bounding box request to OSM's discrete zoom levels.
If you're interested in displaying the OSM tiles over GE, check out Markus Bader's OpenStreetMapLayer.kmz

I've instead displayed the vector version of the trails (linestrings and icons) using a Python script that extracts data from my local OSM Postgresql database. The script runs as a CGI and exports a bounding box worth of data to Google Earth using the KML NetworkLink feature.

Hacky code sample follows (builds the KML Placemarks for the OSM tourism points of interest):

def poisFromBBOX(east,west,south,north):

        result = ""

        conn = psycopg2.connect("dbname=...")

        cur = conn.cursor()

        stmt = cur.mogrify('select name,tourism, st_astext(transform(way, 4326)) from planet_osm_point where tourism is not null and st_contains(st_transform(ST_SetSRID(ST_MakeBox2D(ST_Point(%s, %s), ST_Point(%s, %s)), 4326), 900913), way)', (east, south, west, north))

        cur.execute(stmt)

        tl = cur.fetchall()

        for item in tl:

                name = item[0]               

                if name == None:                        

                    continue

                tourism = item[1]                

                point = item[2]                

                m = re.search('(?<=POINT\()([0-9.]+) ([0-9.]+)', str(point))

                result = result + '<Placemark><styleUrl>#myStyleMap</styleUrl><name>' + escape(str(name))  + '<description>OSM tourism: ' + escape(str(tourism)) + '</description><Point><coordinates>' + str(m.group(1)) + ',' + str(m.group(2)) + ',0</coordinates></Point></Placemark>'

        cur.close()

        conn.close()

        return result