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.

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

Produire un effet de “couleur dégradée” avec Mapserver

Nous avons  travaillé sur un projet de Web mapping dernièrement dans lequel on voulait produire un effet de "couleur dégradée" (fade to white) entre la terre ferme et l'océan. J'ai trouvé dernièrement une superbe carte dans l'"Atlas of Design" (que je recommande fortement pour les maniaques comme moi) de la NACIS (North American Cartographic Information Society).

À la page 64 on trouvera le genre d'effet de "couleur dégradée" en question:

On a l'habitude de voir ce genre d'effet dans un atlas et il s'agit la plupart du temps d'une retouche "photoshop". Il est par contre possible de reproduire l'effet avec une suite de buffers "concentrique" ou "excentrique" autour des entités d'une classe de provinces ou de continents par exemple pour "dégrader" la couleur de l'eau. Je me suis inspiré de d'une recette sur le forum SIG pour le faire dans PostGIS.

La méthode pourrait certe être plus optimale et performante si elle était gérée nativement dans Mapserver (nous acceptons volontier le financement pour le faire), par contre le résultat est superbe.

Donc première étape: chargement... shp2pgsql -s 3163 -c -g the_geom -I -W "latin1" QUARTIER_POLYGON.shp nc_noumea | psql -d nc -h 192.168.6.20 -p 5432 -U nc

Pour créer les anneaux "excentriques" j'ai utilisé ces suites de requêtes SQL. On devra jouer avec la largeur du buffer si on doit afficher le "fade effect" à plusieurs échelles et créer plus d'anneaux. Dans l'exemple suivant (de 0 a 35 mètres) on comprendra que l'effet ne sera visible qu'à très petite échelle. Si on veut appliquer la recette à plus grande échelle, il suffira simplement d'ajouter des anneaux supplémentaires ET à des intervalles plus larges.

NOTE: Merci à Vincent Picavet de la firme Oslandia pour le conseil de l'incroyable commande generate_series(), que je connaissais pas! Vraiment bien cette commande. Cette instruction SQL très compacte, va même créer la table pour moi! drop table nc_noumea_step; create table nc_noumea_step as select 0 as gid,0 as step,the_geom from nc_noumea; insert into nc_noumea_step (gid,step,the_geom) select 1 as gid,5 as step,ST_Difference(ST_Buffer(st_multi,5),ST_Buffer(the_geom,0)) from nc_noumea; insert into nc_noumea_step (gid,step,the_geom) select 2 as gid,10 as step,ST_Difference(ST_Buffer(st_multi,10),ST_Buffer(the_geom,5)) from nc_noumea; insert into nc_noumea_step (gid,step,the_geom) select 3 as gid,15 as step,ST_Difference(ST_Buffer(st_multi,15),ST_Buffer(the_geom,10)) from nc_noumea; insert into nc_noumea_step (gid,step,the_geom) select 4 as gid,20 as step,ST_Difference(ST_Buffer(st_multi,20),ST_Buffer(the_geom,15)) from nc_noumea; insert into nc_noumea_step (gid,step,the_geom) select 5 as gid,25 as step,ST_Difference(ST_Buffer(st_multi,25),ST_Buffer(the_geom,20)) from nc_noumea; insert into nc_noumea_step (gid,step,the_geom) select 6 as gid,30 as step,ST_Difference(ST_Buffer(st_multi,30),ST_Buffer(the_geom,25)) from nc_noumea; insert into nc_noumea_step (gid,step,the_geom) select 7 as gid,35 as step,ST_Difference(ST_Buffer(st_multi,35),ST_Buffer(the_geom,30)) from nc_noumea;
SELECT i as step, ST_Difference(ST_Buffer(st_multi,i),ST_Buffer(st_multi,i-5))
INTO nc_noumea_step
FROM nc_noumea n, generate_series(5,35,5) AS i;

On pourra utiliser un buffer négatif pour produire le même effet mais en dégradant la couleur vers l'intérieur (anneaux concentriques):

...
insert into communes_50_step (gid,cod_commun,step,the_geom) select gid,cod_commun,200 as step,ST_Difference(ST_Buffer(the_geom,-100),ST_Buffer(the_geom,-200)) as the_geom from communes_50;
insert into communes_50_step (gid,cod_commun,step,the_geom) select gid,cod_commun,300 as step,ST_Difference(ST_Buffer(the_geom,-200),ST_Buffer(the_geom,-300)) as the_geom from communes_50;
...

J’ai fait plusieurs essais dont des buffers successifs qui embarquent les uns sur les autres, mais l'effet est moins propre et le résultat moins prévisible. Le meilleur résultat et le plus flexible est vraiment une suite d'anneaux (buffers polygone) successifs qui n'empiètent pas les uns sur les autres. De cette façon on peut gérer la transparence avec couleur de fond pour chaque classe ET/OU des couleurs uniques ET/OU une orthophoto en fond de carte, tout est possible.

Enfin, dans le mapfile j'ai cartographié dans plusieurs CLASS chacun des anneaux désirés avec une transparence spécifique en dégradé :

LAYER
    TYPE POLYGON
    STATUS ON
    GROUP "zone"
    NAME "nc_contour_2_16"
    PROJECTION
        "init=epsg:3163"
    END
    MINSCALEDENOM 4096
    MAXSCALEDENOM 8192
    DATA "DATA_VDN/nc_noumea_step.shp"
    CLASSITEM "step"
    CLASS
        EXPRESSION "0"
        STYLE
            COLOR 255 255 255
        END
    END
    CLASS
        EXPRESSION "5"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 80
        END
    END
    CLASS
        EXPRESSION "10"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 60
        END
    END
    CLASS
        EXPRESSION "15"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 50
        END
    END
    CLASS
        EXPRESSION "20"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 40
        END
    END
    CLASS
        EXPRESSION "25"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 30
        END
    END
    CLASS
        EXPRESSION "30"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 20
        END
    END
    CLASS
        EXPRESSION "35"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 10
        END
    END
END
LAYER
    TYPE POLYGON
    STATUS ON
    GROUP "zone"
    NAME "nc_contour_1_17"
    PROJECTION
        "init=epsg:3163"
    END
    MINSCALEDENOM 2048
    MAXSCALEDENOM 4096
    DATA "DATA_VDN/nc_noumea_step.shp"
    CLASSITEM "step"
    CLASS
        EXPRESSION "0"
        STYLE
            COLOR 255 255 255
            OUTLINECOLOR "#3fc1b7"
            WIDTH 0.4
        END
    END
    CLASS
        EXPRESSION "5"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 50
        END
    END
    CLASS
        EXPRESSION "10"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 30
        END
    END
    CLASS
        EXPRESSION "15"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 10
        END
    END
END
LAYER
    TYPE POLYGON
    STATUS ON
    GROUP "zone"
    NAME "nc_contour_1_18"
    PROJECTION
        "init=epsg:3163"
    END
    MINSCALEDENOM 1500
    MAXSCALEDENOM 2048
    DATA "DATA_VDN/nc_noumea_step.shp"
    CLASSITEM "step"
    CLASS
        EXPRESSION "0"
        STYLE
            COLOR 255 255 255
            OUTLINECOLOR "#3fc1b7"
            WIDTH 0.4
        END
    END
    CLASS
        EXPRESSION "5"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 50
        END
    END
    CLASS
        EXPRESSION "10"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 20
        END
    END
END
LAYER
    TYPE POLYGON
    STATUS ON
    GROUP "zone"
    NAME "nc_contour_1_18"
    PROJECTION
        "init=epsg:3163"
    END
    MINSCALEDENOM 0
    MAXSCALEDENOM 1500
    DATA "DATA_VDN/nc_noumea_step.shp"
    CLASSITEM "step"
    CLASS
        EXPRESSION "0"
        STYLE
            COLOR 255 255 255
            OUTLINECOLOR "#3fc1b7"
            WIDTH 0.4
        END
    END
    CLASS
        EXPRESSION "5"
        STYLE
            COLOR "#3fc1b7"
            OPACITY 50
        END
    END
END

Encodage de base de données Postgresql

Je suis tombé sur un problème d'encodage de ma base de données OSM dernièrement.  Mon besoin était pourtant très simple: je voulais imprimer sur ma carte les noms de rue en MAJUSCULES.  Je pensais donc utiliser la fonction upper() de postgresql via un mapfile, mais étrangement les caractères accentués ne suivaient pas:

osm= select upper('échap');
upper
-------
éCHAP
(1 row)

L'encodage da la base de données était bien UTF-8 puisque tous les caractères accentués apparaissaient convenablement.  Par contre le "Ctype" était à la valeur C. Dans la documentation de Postgresql on y mentionne que:

initdb initializes the database cluster's default locale and character set encoding. The character set encoding, collation order (LC_COLLATE) and character set classes (LC_CTYPE, e.g. upper, lower, digit) can be set separately for a database when it is created. initdb determines those settings for the template1 database, which will serve as the default for all other databases.

postges=l
Name    | Owner    | Encoding  | Collate |    Ctype
osm     | osm      | UTF8      | C       |    C

On a pas vraiment l'habitude de se soucier de cette valeur, mais elle affecte le fonctionnement de certaines fonctions qui traitent les caractères.  Pour contourner le problème, il y a deux options:

  1. Créer une nouvelle COLLATION dans Postgresql;
  2. Recréer le cluster de base de données;

Option 1

Une solution temporaire de contournement serait de créer une nouvelle COLLATION dans la base de données.  Avant on doit s'assurer d'avoir la Collation sur le serveur:

sudo locale-gen fr_CA.UTF-8 en_CA.UTF-8

Ensuite, on peux ajouter une nouvelle Collation correspondant dans la Base de données:

osm= CREATE COLLATION fr_ca_c (LOCALE='fr_CA.utf8');
osm= select upper('échap' COLLATE "fr_ca_c");;
upper
-------
ÉCHAP
(1 row)

Option 2

Cette option plus drastique sera plus robuste et efficiente à long terme.  Elle consiste à recréer le cluster de base de données pour supporter le bon encodage et la bonne Collation(Ctype). En principe, on doit être en mesure de créer la base de données avec les bons réglages de la façon suivante:

sudo su postgres -c'createdb -E UTF8 --lc-ctype en_CA.UTF-8 -T template0 template_postgis'
sudo su postgres -c'createlang -d template_postgis plpgsql;'
sudo su postgres -c'psql -U postgres -d template_postgis -c"CREATE EXTENSION postgis;"'
# # sudo su postgres -c'psql -U postgres -d template_postgis -c"select postgis_lib_version();"'
sudo su postgres -c'psql -U postgres -d template_postgis -c "GRANT ALL ON geometry_columns TO PUBLIC;"'
sudo su postgres -c'psql -U postgres -d template_postgis -c "GRANT ALL ON spatial_ref_sys TO PUBLIC;"'
sudo su postgres -c'psql -U postgres -d template_postgis -c "GRANT ALL ON geography_columns TO PUBLIC;"'
Par la suite, on doit refaire une nouvelle base de données avec ce template et le problème sera réglé: La création d'une nouvelle base de données avec le template UTF-8:
 sudo su postgres
 createdb -E utf8 -T template_postgis osm
 psql -d osm
 postgres= CREATE USER osm WITH PASSWORD 'osm';
 postgres= GRANT ALL PRIVILEGES ON DATABASE osm to osm;
 postgres=\l
 ----------------------------------
 Name | Owner | Encoding | Collate | Ctype
 osm | osm | UTF8 | C | en_CA.UTF-8
 postgres=\q
Il est possible d'avoir des problèmes de création de cette bd J'ai eu des problèmes de création d'une base de données template Postgresql supportant les caractères accentués et/ou l'encodage UTF-8.
sudo su postgres -c'createdb -E UTF8 --lc-ctype en_CA.UTF-8 -T template0 template_postgis'
createdb: database creation failed: ERROR:  invalid locale name en_CA.UTF-8

Comme mentionné plus haut, Postgresql va créer un moteur cluster (service, fichiers config, fichiers de données, logfile, etc) de base de données et ce moteur utilise les configurations d'encodages au moment de l'installation.  Si par mégarde, l'encodage UTF-8 n'était pas disponible, il faut détruire le cluster Postgresql avec pg_dropcluster, mettre à jour les variables locales d'encodage, et en créer un nouveau cluster avec pg_createcluster:

NOTE: Attention vous allez perdre toutes les BD de votre serveur!
sudo /etc/init.d/postgres stop
sudo su postgres
export LANGUAGE=en_CA.UTF-8
export LANG=en_CA.UTF-8
export LC_ALL=en_CA.UTF-8
locale-gen en_CA.UTF-8
pg_dropcluster --stop 9.1 main
pg_createcluster 9.1 main
exit
sudo /etc/init.d/postgres start
Juste pour être certain que tout est ok avec l'encodage de votre user :
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

Chargement Shapefiles dans PostGIS en batch

Il existe un moyen vraiment simple de charger le contenu d'un répertoire de shapefiles dans une base de données Posgresql/PostGIS avec les outils, shp2pgsql et psql. L'utilitaire permettant de charger les commandes SQL (psql) interdit la saisie d'un mot de passe en clair dans une ligne de commande. Pour contourner ce problème, on peut configurer un fichier .pgpass. Le contenu de ce fichier prend cette forme: nom_hote:port:database:nomutilisateur:motdepasse Il est impératif de changer les droits d'accès à ce fichier sinon psql n'en tiendra pas compte.
sudo vim ~/.pgpass
sudo chmod 0600 ~/.pgpass
Par la suite, j'utilise ce petit bash pour générer les lignes de commandes. Ici je transforme le nom du shapefiles en minuscule parce que Postrgresql est "case sensitive" et c'est plus propre de garder les noms de tables en minuscules:
vim shp2sqlfile.sh
---> ajouter ces lignes:
#!/bin/bash
for f in *shp
do
  name=$(basename $f .shp| tr '[A-Z]' '[a-z]')
  echo /usr/lib/postgresql/9.1/bin/shp2pgsql -c -I -W "latin1" -s "EPSG:4326" -g "the_geom" $f $name "| psql -d my_geodb -h localhost -U postgres -w"
done
Ne reste plus qu'a exécuter le batch et rediriger dans un fichier de commandes:
$ sh shp2sqlfile.sh >batch.sh

Créer une zone CANUTEC pour la sécurité civile

Tous ceux qui travaillent en sécurité civile connaissent bien le concepte de la zone CANUTEC.  Cet acronyme est en fait le Centre canadien d’urgence transport qui relève de la Direction générale du transport des marchandises dangereuses (TMD) de Transports Canada.   La zone CANUTEC est simplement une zone de protection établie par les autorités en cas de déversement de matières dangereuses ou d'émanation de gaz toxiques sur le territoire.  Cette zone est établie en fonction de la dangerosité du déversement et de la direction des vents.  On trouve toute la documentation théorique sur le site de transport Canada

Je voulais pouvoir créer une fonction qui allait être appelée depuis différents interfaces.  J'ai donc créé une fonction sur un serveur Posgresql/PostGIS.

CREATE OR REPLACE FUNCTION buildPolyCanutec (in_y DOUBLE PRECISION,
                                             in_x DOUBLE PRECISION,
                                              in_epsg INTEGER,
                                              in_direction DOUBLE PRECISION,
                                              in_dist_isolation integer,
                                              in_dist_aval integer,
                                              in_usr varchar(50),
                                              in_publish boolean)
 RETURNS geometry AS
$BODY$
--/******************************************************************************************
--D Description: Construire un poygone CANUTEC en fonction de point x,y et direction des vents
--A Arguments : in_y = coordonnée x
--              in_x = coordonnée y
--              in_epsg = SRID des coordonnées d'entrées
--              in_direction = direction des vents
--              in_dist_isolation = rayon du buffer circulaire autour du point en mètres
--              in_dist_aval = rayon du buffer circulaire pour montrer la distance à faire en direction des vents
--              in_usr = utilisateur Web qui call la construiction du polygone CANUTEC
--              in_publish = publier la geométrie dans la service Web
--S Sortie : Une géométrie représentant un polygone CANUTEC.
--           Le polygone CANUTEC résulte de l'union de 3 géométries.
--           1)un cercle de diamètre in_dist_isolation
--           2) un cone
--           3) un cercle de diamètre in_dist_aval
-- Voir spec ici http://wwwapps.tc.gc.ca/saf-sec-sur/3/erg-gmu/gmu/modedemploistdiiap.aspx#2
--H Historique: Simon Mercier
--******************************************************************************************/
DECLARE
 geom_buff_isolation geometry ;
 geom_cone_canutec geometry ;
 geom_buff_en_aval geometry ;
 geom_canutec geometry;
 geom_poly_canutec geometry;
 geom_line_translate geometry;
 str_poly_canu varchar(1000);
 x_pts float;
 y_pts float;
 proj_trav integer;
BEGIN

 --On va utiliser le SRID Lambert qui est une projection conforme requise pour concerver les surfaces intactes.
 proj_trav:=32198;

 --On va devoir transformer en mètre les coordonnées.
 select into x_pts ST_XMin(ST_Transform(ST_GeomFromText('POINT('||in_x||' '||in_y||')',in_epsg),proj_trav));
 select into y_pts ST_YMin(ST_Transform(ST_GeomFromText('POINT('||in_x||' '||in_y||')',in_epsg),proj_trav));

 --On va construire un string WKT de la section Cone:
 str_poly_canu:='MULTIPOLYGON (((' || x_pts ||' '|| y_pts+in_dist_isolation ||','|| x_pts + in_dist_aval||' '|| y_pts + (in_dist_aval/2) ||','|| x_pts + in_dist_aval||' '|| y_pts - (in_dist_aval/2) ||','|| x_pts||' '|| y_pts-in_dist_isolation ||','|| x_pts||' '|| y_pts+in_dist_isolation || ')))';

 --On va créer la geom du cone avec le wkt
 select into geom_cone_canutec ST_GeomFromText(str_poly_canu,proj_trav);

 --On a besoin d une geom de l isolation
 select into geom_buff_isolation ST_Buffer(ST_GeomFromText('POINT('||x_pts||' '||y_pts||')',proj_trav), in_dist_isolation);
--On a besoin d une geom en aval
 select into geom_buff_en_aval ST_Buffer(ST_GeomFromText('POINT('||x_pts||' '||y_pts||')',proj_trav), in_dist_aval);

 --On va trimer le cone avec le buffer en aval
 select into geom_canutec ST_Intersection(geom_cone_canutec,geom_buff_en_aval);

 --On va ajouter la zone d isolation
 select into geom_canutec ST_Union(geom_canutec,geom_buff_isolation);

 -- On doit le faire pivoter en fonction des vents... La fonction ST_Rotate de PostGIS utilise un origine
 -- 0,0 pour travailler.  Pour contourner ce problème, j ai trouvé un post sur le net qui fonctionne bien. 
 -- Il faut simplement faire pivoter la géométrie et la déplacer en fonction de son DeltaX et DeltaY. 
 -- On va en premier calculer le DeltaX et DeltaY averc un line string de longueur in_dist_aval
 -- http://geospatial.nomad-labs.com/2007/02/24/rotation-in-postgis/
 select into geom_line_translate ST_GeomFromText('LINESTRING('||x_pts||' '||y_pts||','||x_pts+in_dist_aval||' '||y_pts||')',proj_trav);
 select into geom_line_translate translate( rotate( translate( geom_line_translate, -x(centroid(geom_line_translate)), -y(centroid(geom_line_translate)) ), radians(in_direction)), x(centroid(geom_line_translate)), y(centroid(geom_line_translate)) );

 --On va faire pivoter la géométrie sur son centroïde
 select into geom_poly_canutec translate( rotate( translate( geom_canutec, -x(centroid(geom_canutec)), -y(centroid(geom_canutec)) ), radians(in_direction)), x(centroid(geom_canutec)), y(centroid(geom_canutec)) );
 
 --Et enfin lui appliquer la translation à partir de geom_line_translate.
select into geom_poly_canutec ST_Multi(translate(geom_poly_canutec,x_pts - ST_XMin(ST_StartPoint(geom_line_translate)), y_pts - ST_YMin(ST_StartPoint(geom_line_translate))));
 
 return geom_poly_canutec;
END;
$BODY$
 LANGUAGE 'plpgsql' VOLATILE
;
On va par la suite créer un trigger sur une table de type Linestring pour ajouter automatiquement un polygone CANUTEC dans une table de type Polygon. L'avantage de passer par une table de type Linestring est qu'on pourra facilement extraire une direction et un point de départ à partir de notre ligne.  Je vais utiliser trois fonctions PostGIS pour finaliser le processus: ST_Azimuth(), ST_StartPoint(), ST_EndPoint()

create table p_line_canutec (oid serial, user varchar(100));
select AddGeometryColumn ('public','p_line_canutec','the_geom',32198,'LINESTRING',2);

create table p_poly_canutec (oid serial, user varchar(100));
select AddGeometryColumn ('public','p_poly_canutec','the_geom',32198,'MULTIPOLYGON',2);

Ensuite je vais créer le trigger qui ajoutera automatiquement un polygone CANUTEC dans la table p_poly_canutec à partir de la ligne ajoutée dans la table p_line_canutec. Dans mon exemple j'ai utilisé une zone d'isolation de 8Km et une distance en aval des vents de 50km.

CREATE OR REPLACE FUNCTION auto_build_canutec() RETURNS trigger AS $BODY$
DECLARE
geom geometry;
BEGIN

SELECT INTO geom buildpolycanutec( ST_Y(ST_StartPoint(NEW.the_geom)),
                                   ST_X(ST_EndPoint(NEW.the_geom)),
                                   ST_SRID(NEW.the_geom),
                                   ST_Azimuth(ST_StartPoint(NEW.the_geom),ST_EndPoint(NEW.the_geom))/(2*pi())*360,
                                   8000,50000,NEW._user,'true');

INSERT INTO p_poly_canutec(the_geom,_user) VALUES (geom, NEW._user);

RETURN NEW;
END;
$BODY$
LANGUAGE 'plpgsql';

CREATE TRIGGER trg_build_canutec AFTER INSERT OR UPDATE ON p_line_canutec
FOR EACH ROW EXECUTE PROCEDURE auto_build_canutec();

Dans QGis l'ajout d'une ligne ajoute automatiquement un polygone dans le layer p_poly_canutec.

Greenhouse Gas Emissions Hexagbin Heatmap from MapBox to CartoDB

Version Française

In my last post (Données ouvertes et les émissions de Gaz à effet de serre au Canada) I explained how I convert a CSV data file (from opendata.gc.ca) on greenhouse gases emissions in Canada to build a map. So I produced a tiled map of this project with TileMill deployed on my own server with Wax API.  You can try it now !

I tried to show two things 1) the density of facilities producing those emissions 2) the distribution of GHG emissions in Canada. of course I could used a "Heatmap" to show the density of those value.  But, I tried to produce this density map using a hexagonal grid.

The idea of ​​presenting the amount of GHGs produced in Canada in a points clustered grid came wit this example (Point Clustering) made ​​by the CartoDB developers. More recently, @elcep also proposed a similar approach in his blog. This concept was perfect to illustrate the distribution of greenhouse gases in a 300km x 300km virtual grid.  I builded my layer with this SQL in PostGIS:

SELECT max(install_ges_p.gid) AS gid, sum(install_ges_p.emis_ges) AS gazemis,     st_snaptogrid(install_ges_p.the_geom, 0::double precision, 0::double precision, 300000::double precision, 300000::double precision) AS the_geom_gr FROM install_ges_p GROUP BY st_snaptogrid(install_ges_p.the_geom, 0::double precision, 0::double precision, 300000::double precision, 300000::double precision); I loaded the data from my local "PostgreSQL / PostGIS" into my CartoDB account.  This tool is cool because it's faster to deploy and doesn't require a tiled map production as everything is dynamic.  All we must do is load the table and carto style into your CartoDB account.
@opac_grid:0.1;
@color_gris: #4F1891;
@opac_marker:0.3;
#gridemisgestot[gazemi_nb>0][zoom>=3]{
 marker-fill:@color_gris;
 marker-opacity:@opac_grid;
 marker-line-color:@color_gris;
 marker-line-opacity:@opac_marker;
 marker-allow-overlap:true;
 text-name:"[gazemis]";
 text-face-name:"DejaVu Sans Book";
 text-allow-overlap:true;
 text-halo-radius:1.5;
 text-halo-fill:rgba(255,255,255,0.50);
 text-placement:point;
 [zoom=3]{marker-width:7.5;text-size:4;}
 [zoom=4]{marker-width:15;text-size:9;}
 [zoom=5]{marker-width:30;text-size:18;}
 [zoom=6]{marker-width:60;text-size:36;}
 [zoom=7]{marker-width:120;text-size:36;}
}
Then you can 1) use this API

2) or use the small 'embed_map' provided by CartoDB layer by layer.

https://smercier.cartodb.com/tables/hex_gazemis_cnt/embed_map?sql=SELECT * FROM hex_gazemis_cnt' width='100%' height='480'

Données ouvertes et les émissions de Gaz à effet de serre au Canada

#cop17 ça vous dit quelque chose?  On dénombre pas moins de 60 000 tweets contenant cette mention depuis une trentaine de jours  sur twitter!  Je ne ferai pas ici de la petite propagande climatologique rassurez-vous, par contre c'était bien assez pour piquer ma curiosité ! Trouvant le sujet de cette conférence sur le changement climatique plutôt intéressant d'un point de vue professionnel et personnel, j'ai donc fouillé un peu le sujet pour visualiser la production de  Gaz à effet de serre (GES) au Canada.

En 2008 selon "The World Ressource Institute", le Canada se situait au 8e rang dans l'échelle des plus grands émetteurs de GES avec 561 MtCO2e soit 1.87 % de la production mondiale totale (toutes sources GES confondues).  Globalement c'est peu mais certains diront que c'est déjà trop alors pour voir plus précisément comment se distribuait la production de dioxyde de carbone au Canada j'ai fait quelques recherches sur le Portail de données ouvertes du gouvernement du Canada.  Ce site regorge de données et il suffit de se donner la peine de chercher un peu pour trouver tout plein de données à référence spatiale.  J'ai trouvé des données et une application de cartographie interactive assez rapidement en recherchant le mot clé "Gaz".

Avec ces données téléchargées (format CSV) sur le site d'Environnement Canada, ce que je voulais faire était assez simple: localiser un peu mieux sur une carte la production des GES au Canada.  J'ai géoréférencé les 578 installations émettrices de GES recensées au Canada dans QGis et j'ai ensuite créé des buffers autour des installations pour visualiser leur importance au Canada et un aussi produit grid qui pourra faciliter la représentation cartographique de ces installations.  Parce que je ne voulais faire autre chose qu'une "Heat map", après quelques tentatives,  j'arrive à quelque chose d'intéressant visuellement à l'échelle du Canada (point rouge étant une installation industrielle émettrice de GES et buffer de 100, 250 et 500 km de rayon).

 

Pour produire les valeurs permettant de créer une représentation thématique du sujet (et surtout me faciliter la vie dans mon cas) je transfère les installations et mon grid hexagonal dans Postgresql/Postgis.  J'ai fait quelques requêtes pour comptabiliser 'spatialement' les données.

1) petit nettoyage des données
update hex_emis_ges set "Emissions"='0' where "Emissions"='-99.00000000' or "Emissions" is null
update hex_emis_ges set "Emission_2"='0' where "Emission_2"='-99.00000000' or "Emission_2" is null
update hex_emis_ges set "Emissio3"='0' where "Emissio3"='-99.00000000' or "Emissio3" is null
update hex_emis_ges set emis_ge3_3=to_number("Emissions",'999999') + to_number("Emission_2",'999999')+to_number("Emissio3",'999999');
2) J'ai ensuite simplement représenté l'importance du nombre d'installations émettrices en GES en créant des tables buffers de 100, 250 et 500 km autour des installations.
create table buf500km_install_ges as select "Emissions" as emissions,ST_Buffer(install_ges_p.the_geom,500000,'quad_segs=8') as the_geom from install_ges_p;
CREATE INDEX idx_buf500km_install_ges_geom_gist ON buf500km_install_ges USING gist(the_geom);
create table buf250km_install_ges as select "Emissions" as emissions,ST_Buffer(install_ges_p.the_geom,250000,'quad_segs=8') as the_geom from install_ges_p;
CREATE INDEX idx_buf250km_install_ges_geom_gist ON buf250km_install_ges USING gist(the_geom);
create table buf100km_install_ges as select "Emissions" as emissions,ST_Buffer(install_ges_p.the_geom,100000,'quad_segs=8') as the_geom from install_ges_p;
CREATE INDEX idx_buf100km_install_ges_geom_gist ON buf100km_install_ges USING gist(the_geom);
3) Enfin, j'ai comptabilisé le nombre d'intersections des buffers sur mon grid hexagonal
CREATE TABLE hex_emis_ges_cnt_install_tt as
      select hex_emis_ges.gid, "Latitude", "Longitude" ,"Emissions","Emission_2","Emissio3",
             (select count(*) from buf500km_install_ges where ST_Intersects (buf500km_install_ges.the_geom,
                      hex_emis_ges.the_geom) = true) as nb_buf500,
             (select count(*) from buf250km_install_ges where ST_Intersects(buf250km_install_ges.the_geom,
                      hex_emis_ges.the_geom)=true) as nb_buf250,
             (select count(*) from buf100km_install_ges where ST_Intersects(buf100km_install_ges.the_geom,
                      hex_emis_ges.the_geom)=true) as nb_buf100,the_geom from hex_emis_ges;
CREATE TABLE hex_emis_ges_cnt_install as
select hex_emis_ges_cnt_install_tt.*, (nb_buf500+nb_buf250)+nb_buf100 as nb_tot_install from hex_emis_ges_cnt_install_tt;
drop table hex_emis_ges_cnt_install_tt;
CREATE INDEX idx_hex_emis_ges_cnt_install_geom_gist ON hex_emis_ges_cnt_install USING gist(the_geom);

Le résultat n'est certes pas une révolution ni parfaite et encore moins scientifique, mais l'effet visuel me semble intéressant dans QGis.  La carte montre deux choses:

  • La quantité d'émission de GES par installations industrielles recensées au Canada
  • La concentration et la répartition spatial des installations industrielles au Canada
  Je n'aimais pas vraiment le résultat alors j'ai passé le projet dans la moulinette TileMill pour produire une carte un peu plus soignée.  Il ne me restera qu'à placer la carte tuilée (MBtiles) dans une page pour présenter l'ensemble du territoire et que j'aimerais bien faire avec MapCache cette fois... Prochain billet donc à suivre!  

Autres liens:

http://cait.wri.org/cait.php?page=background&from=yearly http://www.carbon-biodiversity.net/Interactive/CarbonCalculatorNotes http://maps-cartes.ec.gc.ca/indicators-indicateurs/TableView.aspx?ID=15 http://www.protectedplanet.net/  

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