Cartographier des numéros d’autoroutes multiples dans OpenStreetMap

J'ai eu dernièrement à extraire les numéros d'autoroutes des segments de géométries du réseau routier de OpenStreetMap. Il arrive régulièrement que des segments aient plus d'un numéro d'autoroute. Dans ce cas, les numéros seront listés dans le champ 'ref' de la géométrie. Par exemple:
SELECT type, oneway, CHAR_LENGTH(ref) AS reflen,ref
     FROM osm_roads
     WHERE type IN ('motorway', 'trunk')
       AND (ref IS NOT NULL AND ref != '')
AND ref like '%;%'
limit 10;
   type   | oneway | reflen |       ref
----------+--------+--------+-----------------
 trunk    |      1 |      9 | US 412;84
 motorway |      1 |     10 | I 72;US 36
 motorway |      1 |     10 | I 72;US 36
 trunk    |      0 |      9 | NH947;SH6
 motorway |      1 |      9 | I 55;I 69
 trunk    |      1 |     11 | US 10;WI 13
 motorway |      1 |     15 | I 55;I 70;US 40
 trunk    |      0 |     11 | US 278;MS 6
 motorway |      1 |     15 | I 55;I 70;US 40
 motorway |      1 |     11 | I 69;MS 304
(10 rows)
En utilisant un Regex et des array dans la base de données Postgresql/PostGIS, il sera possible d'isoler et multiplier les numéros sur autant de lignes que de numéros trouvés (NOTE: on retrouve une expression régulière semblable dans plusieurs scripts de traitement de données OSM sur le Web): Dans cet exemple ref = 'I 39;I 90'
SELECT 'segment 1' as id,
(regexp_matches('I 39;I 90',
                '([a-z,A-Z]+)?(\s|\-)?([0-9]+)[^;|/|\s]*',
                'g'))[1] as shield_class,
(regexp_matches('I 39;I 90',
                '([a-z,A-Z]+)?(\s|\-)?([0-9]+)[^;|/|\s]*',
                'g'))[3] as shield_no;

    id     | shield_class | shield_no
-----------+--------------+-----------
 segment 1 | I            | 39
 segment 1 | I            | 90
(2 rows)
Dans cet autre exemple ref = 'I 39- 90'
SELECT 'segment 2' as id,
(regexp_matches('I 39- 90',
                '([a-z,A-Z]+)?(\s|\-)?([0-9]+)[^;|/|\s]*',
                'g'))[1] as shield_class,
(regexp_matches('I 39- 90',
                '([a-z,A-Z]+)?(\s|\-)?([0-9]+)[^;|/|\s]*',
                'g'))[3] as shield_no;

    id     | shield_class | shield_no
-----------+--------------+-----------
 segment 2 | I            | 39
 segment 2 |              | 90
(2 rows)
Et enfin ref = 'I 39;US 90/S 109'
SELECT 'segment 3' as id,
(regexp_matches('I 39;US 90/S 109',
               '([a-z,A-Z]+)?(\s|\-)?([0-9]+)[^;|/|\s]*',
               'g'))[1] as shield_class,
(regexp_matches('I 39;US 90/S 109',
                '([a-z,A-Z]+)?(\s|\-)?([0-9]+)[^;|/|\s]*',
                'g'))[3] as shield_no;

    id     | shield_class | shield_no
-----------+--------------+-----------
 segment 3 | I            | 39
 segment 3 | US           | 90
 segment 3 | S            | 109
(3 rows)
En utilisant cette astuce, il sera donc possible de cartographier plusieurs numéros sur un même segment de route

Rechercher des adresses dans Open Street Map

Un étudiant suisse m'a demandé de l'orienter dans sa recherche pour trouver des adresses dans Open Street Map.  Il avait chargé avec succès un fichier de données sources de la Suisse dans sa base de données Postgresql mais avait de la difficulté à trouver les adresses.

On ne trouve pas énormément d'adresses dans OSM car c'est une information très coûteuse à colliger, mais il y en a un peu.   Si vous utilisez l'utilitaire d'importation imposm,  pas de chance, les adresses ne sont pas récupérées par défaut.  Je n'ai jamais essayé, mais il faudrait probablement modifier le script python de paramétrage d'importation defaultmapping.py.   Ne vous en faire pas si vous trouvez la chose un peu compliquée à priori!  Ça l'est réellement et demande du temps pour assimiler l'ensemble des technicalités OSM!

Si vous utilisez l'utilitaire osm2pgsql, vous pouvez récupérer des adresses dans les buildings et les amenities, à condition d'avoir activé l'extension Hstore dans votre base de données, sans oublier de spécifier l'option --hstore dans votre ligne de commande osm2pgsql.  Cette option permet de charger dans la base de données tous les tags et pas seulement ceux spécifiés dans le fichier de style.  Comme j''aime avoir accès à TOUTES la base de données OSM, j'utilise toujours l'option --hstore.  J'ai regardé dans quelques villes en Suisse (dont Lausanne) et effectivement plusieurs bâtiments ont des adresses. Une fois les données chargées avec osm2pgsql et --hstore, voici le genre de requêtes SQL possibles

SELECT osm_id, name, building, shop, amenity,
     tags -> 'tourism'::text AS tourism,
     tags->'denomination' as denomination,
     tags->'religion' as religion,
     tags->'cuisine' as cuisine,
     CAST(tags->'addr:housenumber' as text)||' '||CAST(tags->'addr:street' as text) as address,
     tags->'addr:postcode' as postcode,
     tags->'addr:city' as city,
     tags->'building:levels' as building_levels,
     CASE WHEN ST_IsValid(way) = 'f' THEN           
ST_Buffer(ST_MakeValid(way),0) ELSE way END AS geometry
FROM osm_polygon
WHERE building IS NOT NULL;
Et dans la classe des amenities:
SELECT osm_id, name, amenity, '' as building, shop,
     tags -> 'tourism'::text AS tourism,
     tags->'denomination' as denomination,
     tags->'religion' as religion,
     tags->'cuisine' as cuisine,
     CAST(tags->'addr:housenumber' as text)||' '||CAST(tags->'addr:street' as text) as address,
     tags->'addr:postcode' as postcode,
     tags->'addr:city' as city,
     tags->'building:levels' as building_levels
FROM osm_point
WHERE amenity is not null
OR shop is not null ;
L'autre solution serait d'utiliser l'API Nominatim. Je n'ai pas encore fait d'essais sérieux avec ce service, ici au Canada, ce n'est pas fameux parce que le base de données OSM n'est pas assez complète  Pour terminer il y a aussi le projet OSRM qui pourrait aider à trouver une méthode d'extraction d'adresses.

Un essai avec Vagrant et osmbase

Update le 30 décembre 2013
Suite à un billet de @alanboudrault, sur l'utilisation d'un environnement virtuel pour développer avec Mapserver, j'ai testé l'utilisation de Vagrant pour améliorer notre processus interne de tests unitaires sur nos projets Open Street Map.   J'ai tout de suite aimé la simplicité de l'outil et la rapidité qu'il me donne pour monter un environnement jettable à l'aide de quelques lignes de commandes.  Il faut voir Vagrant comme un outil exploitant un moteur de machine virtuel comme Virtual Box ou VMware Fusion.  Sa force sera d'être rapidement configurable et réutilisable avec un Vagrantfile.

Considérons que vous avez déjà installé Vagrant et VirtualBox. Dans une fenêtre terminale (Mac, Ubuntu ou Windows), on commence par télécharger une box vagrant de type Ubuntu et on va l'initialiser:

mkdir vagrant
cd vagrant
vagrant box add precise32 http://files.vagrantup.com/precise32.box
vagrant init precise32

L'initiation de la VM avec Vagrant, vient de créer un fichier Vagrantfile.  Je garde ce config hyper simple pour l'instant. Il me permettra de me connecter au besoin sur ma VM avec pgAdmin au serveur "localhost" sur le port 5555, et d'accéder au serveur Web localhost via le port 8080 (  http://localhost:8080):

# -*- mode: ruby -*-
# vi: set ft=ruby :

Vagrant.configure("2") do |config|

# Every Vagrant virtual environment requires a box to build off of.
config.vm.box = "precise32"

# Create a forwarded port mapping which allows access to a specific port
# within the machine from a port on the host machine. In the example below,
# accessing "localhost:8080" will access port 80 on the guest machine.
config.vm.network :forwarded_port, guest: 5432, host: 5555
config.vm.network :forwarded_port, guest: 80, host: 8080
...
end
On va démarrer la VM et s'y connecter:
vagrant up
vagrant ssh
Utilisateur Windows: J'utilise ici un client ssh par défaut sur mon macbook. Pour les utilisateurs Windows, je conseille toujours d'utiliser le client ssh MobaXtrem parce qu'il est très pratique pour les moins habitués. L'utilitaire Putty peut aussi être utilisé. Vagrant donnera ce genre d'information (après vagrant up) pour faire la connexion avec MobaXtrem:
Host: 127.0.0.1
Port: 2222
Username: vagrant
Private key: C:/Documents and Settings/mapgears/.vagrant.d/insecure_private_key

J'utilise toujours le user par défaut de vagrant. Un fois connecté au serveur, configurer la langue du nouveau serveur dans le profile vagrant. Ça reste une étape toujours importante pour supporter convenablement les accents dans votre base de données Postgresql:

export LANGUAGE=en_CA.UTF-8
export LANG=en_CA.UTF-8
export LC_ALL=en_CA.UTF-8
sudo localedef -i en_CA -f UTF-8 en_CA
sudo locale-gen en_CA.UTF-8
sudo update-locale
sudo sudo dpkg-reconfigure locales
On doit maintenant faire la mise à jour du serveur avec un script VirtualBox et brancher le dépôt logiciels géomatiques UbuntuGIS (version precise):
sudo sh postinstall.sh
sudo apt-get update
sudo apt-get install python-software-properties
sudo add-apt-repository ppa:ubuntugis/ppa
sudo apt-get update
Installation de Postgresql/PostGIS
sudo apt-get install -y postgresql-9.1 postgresql-server-dev-9.1 postgresql-contrib-9.1 postgis
Installation d'outils
sudo apt-get install -y apache2 binutils gdal-bin checkinstall git vim screen make python-virtualenv python-pip python-all-dev osm2pgsql osmosis 
Installation Mapserver / Mapcache
sudo apt-get install -y cgi-mapserver mapserver-bin libmapcache mapcache-cgi mapcache-tools libapache2-mod-mapcache
sudo mkdir /tmp/ms_tmp
sudo chown www-data:www-data /tmp/ms_tmp
NOTE: Pour faire toutes les opérations en une seule étape, créer un fichier init.sh avec vi à l'aide de ce Gist init.sh

Par la suite on installe un nouvelle BD Postgresql/PostGIS OSM avec un encodage UTF-8 (voir ce billet pour instructions):

Ne reste plus qu'à faire un peu de magie avec ce repo osmbase. Le makefile contenu dans ce projet va télécharger les données requises et lancer imposm avec notre petit fichier OSM préféré de la Nouvelle-Calédonie pour construire ma carte QA et de plus j'active un environnement virtuel python. Enfin, simplement modifier les variables de départ du makefile pour vos besoins.

virtualenv venv
source venv/bin/activate
pip install imposm
git clone https://github.com/smercier/osmbase.git
cd osmbase
make
make imposm
Si l'installation se passe sans erreur, la URL suivante devrait fonctionner! http://localhost:8080/cgi-bin/mapserv?map=/home/vagrant/osmbase/osmbase.map&LAYERS=default&template=openlayers Une fois le QA terminé, on a qu'a détruire la VM et recommencer le processus à nouveau pour un autre besoin.  Très simple, j'aime vraiment :-)

Cartographie d’un réseau routier avec MapServer: savoir gérer l’ordre d’affichage

Lorsqu’on veut cartographier un réseau routier, l’ordre d’affichage des routes est important pour que les voies les plus importantes ne soient pas coupées par des voies secondaires, par exemple.

Image 1 - résultat désiré

Image 2 - résultat brut

L'ordre d'affichage dépend de plusieurs facteurs.  Cet article explique le comportement de MapServer face à l'ordre d'affichage des entités géographiques afin de vous permettre d'obtenir un résultat "propre" dans l'affichage de votre réseau routier.  Il est important de noter que certaines des règles qui suivent ne s'appliquent qu'aux LAYERs de type LINE.

1) Premier lu - Premier dessiné

Le premier élément à considérer est, qu'à l'intérieur d'un même layer, les éléments sont dessinés selon l'ordre de lecture de la donnée indépendamment de la CLASS à laquelle ils appartiennent. En d'autres mots, le premier segment du réseau routier lu sera le premier dessiné. Il est donc nécessaire de classer la source. Par exemple, avec une BD Postgresql/PostGIS, on peut classer les données selon un z-index (z_order dans Open Street Map) et selon la longueur du segment :

DATA "geometry from (select osm_id, geometry, name, type
      from osm_roads
      where type in ('secondary', 'tertiary', 'residential')
      order by z_order asc, st_length(geometry) asc)
      as foo using unique osm_id using srid=900913"

Si vous utilisez un shapefile, l'utilitaire sortshp pourrait être très utile. 

NOTE: Cette règle s'applique à tous les types de layers.

2) Nombre de style dans la classe

Le nombre de STYLE dans une CLASS influence l'ordre d'affichage. Lorsqu'à l'intérieur d'un LAYER, des CLASS comprennent plus qu'un STYLE, le premier STYLE de chaque CLASS est dessiné avant de procéder au 2e STYLE des CLASS et ainsi de suite.

3) Outline en premier

Lorsque les premiers STYLE des CLASS comportent un OUTLINEWIDTH, les outlines de toutes les CLASS sont dessinés avant l'intérieur.  Cet aspet spécifique n'est pas le fruit du hasard et a été implanté dans la version 5.4 et documenté dans la RFC-49.

Exemple: mapfile incorrect

CLASS  # routes secondaires
  EXPRESSION /secondary/
  STYLE
    WIDTH 12
    COLOR 255 0 0 # rouge
    OUTLINEWIDTH 2
    OUTLINECOLOR 153 111 57 # brun
  END
  STYLE
    WIDTH 11
    COLOR 223 197 124 # orange
  END
END

CLASS # routes tertaires
  EXPRESSION /tertiary/
  STYLE
    WIDTH 11
    OUTLINEWIDTH 1
    OUTLINECOLOR 193 188 157 # gris
    COLOR 255 253 139 # jaune
  END
END

CLASS # routes résidentielles
  STYLE
    WIDTH 11
    OUTLINEWIDTH 1
    OUTLINECOLOR 103 181 157  # vert
  END
  STYLE
    WIDTH 11
    COLOR  238 225 226  # rose
  END
END

En considérant que la source a été triée dans cet ordre: résidentielles, tertiaires, secondaires, l'ordre d'affichage se déroule comme suit:

  1. Outline (vert) des routes résidentielles
  2. Outline (gris) des routes tertiaires
  3. Outline (brun) des routes secondaires
  4. Premier style (uniquement un outline) des routes résidentielles
  5. L'unique style (jaune) des routes tertiaires
  6. Premier style (rouge) des routes secondaires
  7. Deuxième style (rose) des routes résidentielles
  8. Deuxième style (orange) des routes secondaires
 

Exemple: mapfile correct

CLASS  # routes secondaires
  EXPRESSION /secondary/
  STYLE
    WIDTH 12
    COLOR 255 0 0 # rouge
    OUTLINEWIDTH 2
    OUTLINECOLOR 153 111 57 # brun
  END
  STYLE
    WIDTH 11
    COLOR 223 197 124 # orange
  END
END

CLASS # routes tertaires
  EXPRESSION /tertiary/
  STYLE
    WIDTH 11
    OUTLINEWIDTH 1
    OUTLINECOLOR 193 188 157 # gris
  END
  STYLE
    WIDTH 11
    COLOR 255 253 139 # jaune
  END
END

CLASS # routes résidentielles
  STYLE
    WIDTH 11
    OUTLINEWIDTH 1
    OUTLINECOLOR 103 181 157  # vert
  END
  STYLE
    WIDTH 11
    COLOR  238 225 226  # rose
  END
END

imposm and PostGIS 2.0

I recently load OSM planet in PostGIS 2.0.  It works but you have to created manually your database:
 sudo su postgres
 createdb -E utf8 -O postgres osmplanet
 psql -d osmplanet -c "CREATE EXTENSION postgis;"
 createuser --no-superuser --no-createrole --createdb osm
 psql -d osmplanet -c "ALTER USER osm WITH PASSWORD 'osm';"
 psql -d osmplanet -f /usr/local/lib64/python2.7/site-packages/imposm/900913.sql
 psql -d osmplanet -c "ALTER TABLE geometry_columns OWNER TO osm;”
 psql -d osmplanet -c "ALTER TABLE spatial_ref_sys OWNER TO osm;"
 sudo vim /var/lib/pgsql/data/pg_hba.conf
 -->add this line "host osmplanet osm  127.0.0.1/32 md5"

 sudo virtualenv venv
 source venv/bin/activate
 venv/bin/pip install imposm
I use multiprocessing option with "imposm" in one line command read/write
 sudo imposm -c 8 --overwrite-cache --read --write --database osmplanet planet-latest.osm.bz2
Final database is 66 GB at this moment
postgres@osm:/home/smercier> psql -c "SELECT pg_size_pretty(pg_tablespace_size('pg_default'));" osm
 pg_size_pretty
----------------
 66 GB
(1 row)

Imposm pour Open Streets Map avec AQréseau dans votre SIG bureautique

J'ai utilisé pour un client, l'utilitaire imposm dans le but de créer une carte de base spécifique pour ses besoins. Ce logiciel Open source s'exécute dans un environnement virtuel python "virtualenv". J'ai été surpris par l'efficacité de l'outil et surtout l'énorme potentiel qu'il offre pour avoir accès à une donnée de contexte gratuite, disponibles en format vectoriel de façon uniformes pour la planète !  Cet utilitaire a le mérite de grandement simplifier l'utilisation des données du monstre OSM en créant un modèle de donnée plus proche du SIG que du GPS (modèle par défaut de OSM). Le découpage des données OSM proposé par impom est documenté sur le site.

Pour utiliser cet outil vous avez besoin:
  • d'un Mac ou Ubuntu
  • d'un serveur Postgresql / PostGIS  ( en mode local sur votre Laptop ou PC pas de problème)
  • des librairies et outil imposm  (l'installation est simple, suivre la procédure sur le site Web)
example de données disponibles

J'ai utilisé un fichier de données OSM découpé selon les limites de la province de Québec à partir du serveur Geofabrik.  Geofabrik est une entreprise d'Allemagne qui offre des services professionnelles aux entreprises et organismes publics désirant utiliser l'information Open Streets Map pour leurs besoins.  Sur ce serveur, les données sont mises à jour de façon quotidienne.  Le format utilisé par IMPOSM est le pbf (Protocolbuffer Binary Format ) qui est plus compressé et plus rapide.  Dans mon cas, j'avais déjà un serveur Postgresql / PostGIS, donc je n'ai eu qu'à créer une nouvelle base de données.  J'ai ajouté dans ma procédure, l'ajoute de la Géobase AQréseau

1)Installer imposm
sudo apt-get install build-essential python-dev protobuf-compiler 
                     libprotobuf-dev libtokyocabinet-dev python-psycopg2 
                     libgeos-c1
sudo apt-get install libgeos-dev
sudo apt-get install python-virtualenv
virtualenv venv
source venv/bin/activate
pip install Shapely==1.2.10
pip install https://github.com/olt/shapely_speedups/tarball/master
pip install imposm
2)Télécharger un fichier de données OSM via Geofabrik
wget http://download.geofabrik.de/osm/north-america/canada/quebec.osm.pbf
3)Créer une base de données PostGIS  'osm'
createdb -E utf8 -O postgres -T postgis osm
4)Charger les données avec l'outil imposm
imposm --read quebec.osm.pbf
imposm --write --database osm --host localhost --user osm
imposm  --optimize -d osm
5)Je vais en profiter pour ajouter la Géobase d'Adresses Québec à mon serveur  de cette façon:
wget http://adressesquebec.gouv.qc.ca/zip/aqreseau.zip
unzip aqreseau.zip
cd ../ESRI(shp)/
shp2pgsql -s 3798 -c -g the_geom -I -W "latin1" AQ_ROUTES.shp aq_routes >aq_routes.sql
export PGCLIENTENCODING=LATIN1
psql -d osm -U osm -W -f aq_routes.sql
NOTE: Attention au système de projection de la base de données epsg:3798

Pour ce qui est de mon importation de la province, la liste des tables du modèle, avec le nombre d'enregistrements pour la province de Québec est listé ici.  De plus, on remarquera que toutes les tables on un préfixe (osm_new).  Cette particularité est liée aux fonctionnalités facilitant les mises à jour des données OSM (mise à jour différentielle non supportée).

osm=# SELECT schemaname,relname,n_live_tup 
FROM pg_stat_user_tables 
ORDER BY n_live_tup DESC;
schemaname |         relname          | n_live_tup
------------+--------------------------+------------
public     | aq_routes                |     408276
public     | osm_new_waterways        |     245797
public     | osm_new_minorroads       |     195209
public     | osm_new_waterareas       |     131509
public     | osm_new_landusages       |      98800
public     | osm_new_mainroads_gen1   |      69004
public     | osm_new_mainroads_gen0   |      69004
public     | osm_new_mainroads        |      69004
public     | osm_new_waterareas_gen1  |      47531
public     | osm_new_buildings        |      38948
public     | osm_new_landusages_gen1  |      34800
public     | osm_new_landusages_gen0  |      11318
public     | osm_new_waterareas_gen0  |       9585
public     | osm_new_motorways_gen0   |       7320
public     | osm_new_motorways_gen1   |       7320
public     | osm_new_motorways        |       7320
public     | osm_new_railways_gen1    |       6129
public     | osm_new_railways_gen0    |       6129
public     | osm_new_railways         |       6129
public     | osm_new_transport_points |       5010
public     | osm_new_places           |       3191
public     | osm_new_amenities        |       1933
public     | osm_new_admin            |        473
public     | osm_new_transport_areas  |        222
public     | osm_new_aeroways         |        193

La taille de la base de données OSM pour la province de Québec n'est que de 635Mb dans la base de données si on soustrait la taille des données Adresses Québec (AQréseau).

osm=# SELECT pg_size_pretty(pg_tablespace_size('pg_default'));
 pg_size_pretty
----------------
 744 MB
(1 row)
osm=# SELECT pg_size_pretty(pg_relation_size('aq_routes')) as table_size;
 table_size
------------
 109 MB
(1 row)

Du coup, je me retrouve avec une base de données de contexte pratique faciles à utiliser avec QGis ou n'importe quel autre SIG bureautique pouvant se connecter a Postgersql / PostGIS.

Ajoute d'une couche de parc dans mon projet QGis

MapOSMatic

MapOSMatic is a web services to generate maps of cities or towns,  including index of streets, from OpenStreetMap data. From About project, MapOSMatic was launched with an idea of Gilles Lamiral, a contributor to OpenStreetMap and Free Software in the region in Rennes France. From his idea, a group of hackers crazy met Hackfest during one week in August 2009 and has transformed the idea into reality Gilles. The group of hackers crazy wishes to thank Gilles for sharing this brilliant idea! It is made of two components:
  • maposmatic, the web front-end. An application written using the Django framework allows to submit and visualize map rendering jobs. The rendering is done in the background by a daemon called maposmaticd;
  • ocitysmap, the back-end that generates the map. It is available as a Python module, used both by the maposmatic daemon (above) and by a sample command line application.
You are invited to contribute to the project on  Savannah