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

Atelier initiation à OpenStreetMap et QGis

Du 17 au 20 juin dernier avait lieu le 48e colloque annuel de l’Association des cartothèques et archives cartographiques du Canada (ACACC).   Les organisateurs sont entrés en contact avec Mapgears pour organiser un atelier sur OpenStreetMap et QGis.   Le défi était tout de même relativement important pour plusieurs participants qui n'avaient pas encore eu la chance d'ouvrir le logiciel QGis.  Comme promis aux participants, voici donc les documents PDF de l'atelier.

Scribeui 0.7 beta

En 2013, Mapgears lance le projet ScribeUI destiné à simplifier la création de cartes pour le Web basé sur Mapserver et l'intégration du langage scribe pour mapserver.  L'objectif était vraiment simple: améliorer notre efficacité d'écriture d'un mapfile pour produire des cartes multiniveau tout en conservant la grande flexibilité de Mapserver.  La première version du projet était vraiment intéressante et on décide de faire un Google summer of code pour améliorer le design de l'outil.  Jessica Lapointe (de l'équipe Mapgears) avait d'ailleurs eu de bons mots à propos de son projet ScribeUI au GsoC 2013 provenant de la responsable de la fondation OSGeo.

Suite à des travaux importants réalisés chez Mapgears, je teste depuis plusieurs jours, une nouvelle version de ScribeUI et je suis vraiment épaté du résultat. Le projet n'est pas encore sur Github, mais il y sera prochainement.  On a remodelé l'interface, ajouté de nouvelles fonctions mais surtout changé le framework Flask pour Pyramid. L'application est maintenant beaucoup plus stable même si quelques petits bogues persistent. La bonne nouvelle est que Jessica Lapointe poursuivra bientôt les travaux entamés sur l'application beta pour améliorer l'utilisation générale du produit.

Personnellement, je ne peux plus me passer de l'outil pour produire des cartes web.  ScribeUI offre une méthode rapide, structurée et efficace pour faire un projet Web Mapping.  Combiné à Git, je peux transférer un projet d'une VM local à un serveur Cloud en quelques minutes, c'est vraiment pratique.

.

En plus de la syntaxe scribe très pratique, deux nouvelles fonctions vraiment cool arrivent avec cette version!  On a intégré un processus Git à la production carto avec ScribeUI. On se disait que produire une carte web est un peu comme faire de la programmation et que git mériterait d'être plus largement utilisé dans ce contexte.  C'est aussi une bonne pratique pour structurer le travail à faire, suivre l'évolution des versions de la carte et faciliter le transfert du projet cartographique d'un environnement à l'autre (DEV, TEST, PROD).  Directement dans l'application ScribeUI, on pourra maintenant, cloner puller et pusher toutes la config d'un projet carto sur un Git.  J'ai un projet complet (data et mapping ) pouvant être cloné ici : https://github.com/smercier/gitscribe.

On a ajouté un lien direct sur le fichier "readme.md" de votre projet pour favoriser le réflex de documentation. Deux clics suffisent pour éditer le "readme.md" de votre git avec la syntaxe mardown.  Pour ma part, c'est dans le readme que je documente les grandes lignes du makefile d'installation de données (quand c'est requis).

On a amélioré la coloration syntaxique et on pourra mieux apprécier la nouvelle option pour commenter de grande section de mapfile avec un code ouvrant (/*) et fermant  (*/) , ça j'aime VRAIMENT.  On a ajouté un lien direct sur l'aide en ligne(http://mapserver.org/mapfile/)  des mapfiles keywords en cliquant sur ALT et le keywords.

Écrire des mapfiles avec des variables et des codes de niveau de zoom, qui correspond bien plus à notre nouvelle réalité de production cartographique pour OpenLayers ou Leaflet.  Faire une carte avec les données OpenStreetMap devient soudainement bien plus plaisant.

L'autre fonction nouvellement ajoutée est un UI pour lancer des job Mapcache.  Le panorama permet de spécifier le Zoom levels à produire, le Metatile size, le Grid disponible et un Extent spécifique.  Ce dernier peut être dérivé du Map Extent de mapfile mais aussi à partir shapefile ou même d'une requête provenant d'un serveur Postgis.   L'intégration de nouvelles options sera documentée sur github.

J'ai aussi testé l'application avec des Mapfiles standard très volumineux (pour une carte maritime avec des données S-57) 5 mapfiles de 2000 lignes et plus dont un de 3042 lignes.  La navigation entre les groupes me permette de passer rapidement d'un include mapfile très rapidement et rend l'édition de la carte bien plus rapide.

À suivre donc, d'autres infos suivront prochainement.  Les idées ne manquent pas pour ce projet.

Plug-in Oracle pour Mapserver dans UbuntuGIS

Mapgears a réalisé un plug-in Oracle pour Mapserver et GDAL/OGR qui simplifie énormément la configuration d'une configuration Mapserver pour Oracle.

Installer le client Oracle

Oracle met à la disposition des organisations un portail regroupant une panoplie d’outils et de tutoriels logiciels Open Source d’Oracle. Grâce à ce portail, il est possible d’installer une version gratuite de la base de données Oracle et le client Oracle XE plus facilement sur un poste Linux.  Par contre, actuellement les packages de ce type d'installation datent de 2006 et il ne s'agit pas de la dernière version du client.  Ça fonctionne, mais l'option manque un peu de rigueur pour monter une infrastructure solide pour une organisation en mode "stable".

Pour permettre la connectivité entre Mapserver et Oracle, je recommande l'installation de la version du client Oracle officiellement supportée:  Oracle InstantClient.  Vous trouverez sur le site Web d'Oracle toutes les versions OS supportées.  Dans cette procédure, j'utilise les sources (dans fichier zip) et non les packages Redhat (rpm).  Pour télécharger les sources officielles, il faudra vous enregistrer sur le site Web d'Oracle et accepter les termes de la licence.  Par la suite, faire l'installation via la méthode suivante :

  • Télécharger les fichiers zip requis pour l'installation de la version Linux x86 InstantClient v 11.2. Les fichiers zip requis seront les suivants :
  1. Instant Client Package - Basic Lite ( instantclient-basic-linux-11.2.0.4.0.zip )
  2. Instant Client Package - SDK ( instantclient-sdk-linux-11.2.0.4.0.zip )
  3. Instant Client Package - SQL*Plus ( instantclient-sqlplus-linux-11.2.0.4.0.zip )
  • Pour installer le client, aucune installation n'est requise.  Il faut simplement dézipper les fichiers sources dans le répertoire de votre choix (/opt/ dans mon cas):
sudo unzip instantclient-sqlplus-linux-11.2.0.4.0.zip -d /opt/
sudo unzip instantclient-basic-linux-11.2.0.4.0.zip -d /opt/
sudo unzip instantclient-sdk-linux-11.2.0.4.0.zip -d /opt/
ls -l /opt/instantclient_11_2/
-rwxrwxr-x 1 root root     22004 Aug 25 11:18 adrci
-rw-rw-r-- 1 root root       437 Aug 25 11:18 BASIC_README
-rwxrwxr-x 1 root root     38461 Aug 25 11:18 genezi
-r-xr-xr-x 1 root root       368 Aug 25 11:18 glogin.sql
-rwxrwxr-x 1 root root  44316855 Aug 25 11:18 libclntsh.so.11.1
-r-xr-xr-x 1 root root   7098468 Aug 25 11:18 libnnz11.so
-rwxrwxr-x 1 root root   1881900 Aug 25 11:18 libocci.so.11.1
-rwxrwxr-x 1 root root 118729922 Aug 25 11:18 libociei.so
-r-xr-xr-x 1 root root    152585 Aug 25 11:18 libocijdbc11.so
-r-xr-xr-x 1 root root   1499831 Aug 25 11:18 libsqlplusic.so
-r-xr-xr-x 1 root root   1216654 Aug 25 11:18 libsqlplus.so
-r--r--r-- 1 root root   2091135 Aug 25 11:18 ojdbc5.jar
-r--r--r-- 1 root root   2739616 Aug 25 11:18 ojdbc6.jar
drwxrwxr-x 4 root root      4096 Aug 25 11:18 sdk
-r-xr-xr-x 1 root root      6913 Aug 25 11:18 sqlplus
-rw-rw-r-- 1 root root       441 Aug 25 11:18 SQLPLUS_README
-rwxrwxr-x 1 root root    160203 Aug 25 11:18 uidrvci
-rw-rw-r-- 1 root root     66779 Aug 25 11:18 xstreams.jar
  • Pour configurer le client par la méthode couramment utilisée par la majorité des organisations, on va créer l'arborescence de répertoires network et y ajouter un fichier tnsnames fourni par le DBA. ( cette étape est optionnelle)
sudo mkdir /opt/instantclient_11_2/network
sudo mkdir /opt/instantclient_11_2/network/admin
sudo vim /opt/instantclient_11_2/network/admin/tnsnames.ora
  • Pour avoir le runtime library path du le Client Oracle pour tous les utilisateurs, il est préférable d’ajouter un fichier de configuration dans le répertoire de chargement de librairie, pour que tous les logiciels puisent y accéder.
sudo vim /etc/ld.so.conf.d/oracle.conf
/opt/instantclient_11_2 ==>ligne à ajouter dans le fichier
sudo ldconfig
  • Validation de connexion direct (sans TNS)
sqlplus [utilisateur]/[mot_passe]@//[nom_server_ou_ip]:1521/[SID]
  • Validation de connexion (avec TNS)
sqlplus [utilisateur]@[SID]

Installer Mapserver et le plugin Oracle

Les instructions de base pour installer le plugin Oracle pour Mapserver sont maintenant documenté sur le Trac UbuntuGIS.
  • 1_ Connecter les packages privés de UbuntuGIS
sudo apt-get install python-software-properties
sudo add-apt-repository ppa:ubuntugis/ppa
sudo apt-get update
  • 2_ Installer Mapserver
sudo apt-get install cgi-mapserver mapserver-bin libmapcache mapcache-cgi mapcache-tools libapache2-mod-mapcache libmsplugin-oracle-src
  • 3_ Activer le plug-in Oracle pour Mapserver.  NOTE 1: on doit créer des liens symboliques avec les librairies requises par le plug-in.  Oracle ne respecte pas la nomenclature et les façons de faire Debian alors on doit s'adapter un peu.  NOTE 2: la variable d'environnement n'est pas requise pour rouler Mapserver une fois le plug-in activé, elle est seulement requise pour l'activation. NOTE 3: Attention, dans mon cas, je me suis retrouvé avec une librairie 11.1 dans un package annoncé 11.2! Je n’ai pas de félicitation à faire ici...
sudo ln -s /opt/instantclient_11_2/sdk/include/ /opt/instantclient_11_2/include

sudo ln -s /opt/instantclient_11_2/libocci.so.11.1 /opt/instantclient_11_2/libocci.so

sudo ln -s /opt/instantclient_11_2/libclntsh.so.11.1 /opt/instantclient_11_2/libclntsh.so

sudo ORACLE_HOME=/opt/instantclient_11_2 mapserver-oracle-build yes
  • 4_ Par la suite utiliser un répertoire relatif dans le paramètre PLUGIN du LAYER ou paramétriser le répertoire dans l'entête du mapfile avec "CONFIG MS_PLUGIN_DIR"
LAYER
  ...
  CONNECTIONTYPE PLUGIN
  PLUGIN "/usr/lib/msplugins/libmsplugin_oracle.so"
  DATA "GEOM FROM ZZ_EOLIENNES_P USING UNIQUE EOLIENNES_ID SRID 4326"
  ...
END
Note: jeu de données test

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 :-)

GEOTRANSFORM et Mauve Nanane …

En début de semaine, j'ai pris connaissance des travaux de mon collègue @AlanBoudreault sur une nouvelle fonctionnalité vraiment intéressante qui sera ajoutée dans la version Mapserver 6.4.  Le RFC-89  ajoutera de nouvelles fonctions de simplification, qui seront effectuées directement par le moteur cartographique Mapserver sans intermédiaire.

J'en ai profité en même temps pour faire une critique cinglante (à l'interne mais c'était humoristique quand même  :-) ) sur l'horrible Mauve nanane du site Web officiel de mapserver.org. Je sais, tout est question de goût, reste que je pense qu'il y trop de mauve dans la mise en page et qu'on devrait faire disparaître tous ces en-têtes mauves qui me font toujours sursauter un peu.

Ce qui m'a valu un retour tout aussi cinglant de la part de bon ami @dmorissette, qui m'a proposé un nouveau titre pour mon blog! Pour ma défense, c'est un paramètre par défaut de mon thème WordPress, mais promis je vais chercher à solutionner ladite chose.

Une carte choroplèthe trigrid

Après avoir créé la carte sur les gaz à effet de serre, je voulais refaire l'expérience avec un grid en forme de triangle, sachant que l'utilisation d'un grid rectangulaire ou hexagonal donne un superbe effet pour cartographier une choroplèthe. Après quelques recherches, j'ai trouvé cet exemple qui s'approchait à ce que je voulais reproduire:

J'ai donc récupéré quelques données géométriques sur les divisions de recensement et des données statistiques sur le site Web de Statistique Canada pour créer une carte Choroplèthe selon les divisions de recensement canadiennes.

Après quelques manipulations dans QGis, j'arrive rapidement à créer une choroplèthe, tout ce qu'il y a de plus standard et habituelle (pour ne pas dire sans saveur ni odeur) basée sur une régression quantile en 5 classes:

Même que, pour mieux représenter la répartition de la population, je préfère utiliser un écoumène de Statistique Canada.  Ça donne un bien meilleur résultat et montre mieux l'occupation du territoire.

J'ai aussi testé une requête SQL proposée par @tokumin sur stackoverflow permettant de calculer directement dans Postgresql la même régression Quantile. Le résultat basé sur l'utilisation de la clause OVER (une windows function bien expliquée dans ce billet) est vraiment intéressant et donne un résultat presque identique à QGis.

SELECT ntile, avg(nb) AS avgAmount, max(nb) AS maxAmount,
       min(nb) AS  minAmount
FROM (SELECT nb, ntile(5) OVER (ORDER BY nb) AS ntile
FROM cdec_p_4326) x
GROUP BY ntile ORDER BY ntile<code>

-- ntile avgamount maxamount minamount
-- 1 11736.7796610169 18000 629
-- 2 23135.5593220339 31138 18036
-- 3 39888.5762711864 50512 31333
-- 4 72870.1206896552 104075 50900
-- 5 428264.25862069 2615060 105719</code>

Pour reproduire  l'effet désiré, j'ai placé dans une boucle la création de 8 triangles à l'aide de quelques formules mathématiques sur les triangles de cette façon:

Je répète ensuite la création de cet ensemble de triangles sur toute la largeur de l'étendue souhaitée et empile les uns par-dessus les autres, ce même ensemble de triangles, sur toute la hauteur voulue. Vous trouverez les quelques lignes de code SQL permettant de créer le trigrid avec Postgresql/PostGIS dans mon github.

À l'aide de cette suite d'instructions et de la fonction trigrid, on pourra créer deux classes d'entités géométriques permettant de préparer une carte choroplèthe trigrid :
-- créer une table tirgrid level 1
drop table trigrid_l1;
create table trigrid_l1(gid integer,tx integer, ty integer);
SELECT AddGeometryColumn ('','trigrid_l1','the_geom',4326,'POLYGON',2);
select trigrid (-150,36,100,44,0.13,4326,'trigrid_l1');
CREATE INDEX sidx_trigrid_l1 ON trigrid_l1 USING GIST ( the_geom );
CREATE UNIQUE INDEX idx_trigrid_l1 ON trigrid_l1 (gid);
select count(*) from trigrid_l1;
-- 605200

-- créer une table tirgrid level 2
drop table trigrid_l2;
create table trigrid_l2(gid integer,tx integer, ty integer);
SELECT AddGeometryColumn ('','trigrid_l2','the_geom',4326,'POLYGON',2);
select trigrid (-150,36,100,44,0.065,4326,'trigrid_l2');
CREATE INDEX sidx_trigrid_l2 ON trigrid_l2 USING GIST ( the_geom );
CREATE UNIQUE INDEX idx_trigrid_l2 ON trigrid_l2 (gid);
select count(*) from trigrid_l2;
-- 2410968

Enfin, pour simplifier la création de la thématique avec Mapserver, j'ai transféré ma classification quantile des divisions de recensement sur mes deux trigrid (level 1 et 2), par une association avec ce genre de requête:

ALTER table trigrid_l1 add column cls integer;
update trigrid_l1 SET cls=5 FROM cdec_p_4326 WHERE ST_Intersects( trigrid_l1.the_geom,cdec_p_4326.the_geom) and cdec_p_4326.nb >105719;
--
update trigrid_l1 SET cls=4 FROM cdec_p_4326 WHERE ST_Intersects( trigrid_l1.the_geom,cdec_p_4326.the_geom) and cdec_p_4326.nb >50900 and cdec_p_4326.nb ...
Au final, j'aime bien et ça donne un résultat pas mal intéressant! Essayez-le http://smercier.github.com/trigrid/trigrid.html!

619 fonctions avancées d’analyses spatiales

QGis n'est pas un outil très utilisé au Québec, par contre les Analystes en Géomatique doivent absolument y jeter un oeil.  Nous utilisons QGis dans la très grande majorité de nos projets chez Solutions Mapgears.  Je n'avais pas encore regardé ce que pouvait offrir l'extension Sextante (aussi disponible sous ArcGIS et gvSIG), probablement parce que ma todolist est un peu trop longue! Bref, ce post d'Oslandia n'avait pas eu encore de suite.

Incroyable, mais on retrouve dans la toolbox de l'extension, pas moins de 619 géoalgorithmes!  GDAL, GRASS, SAGA, OTB(télédétection), R (statistique)....     Vrai que l'interface utilisateur de GRASS n'est pas très sexy, mais on pourra quand découvrir sa puissance via la liste interminable de fonctions dans le toolbar Sextante.  De plus, l'outil dispose d'un "Modeler" mais là vraiment je peux pas en dire plus je ne l'ai pas testé ...

J'ai aussi découvert l'outil d'analyses géoscientifiques System for Automated Geoscientific Analyses (SAGA) un logiciel FOSS.  Le plugin Sextante propose 279 algo de cette plateforme. Le développement semble avoir cessé en 2007, par contre l'outil est le fruit de plusieurs années de développement et vaux vraiment la peine d'être examiné. À mettre sur votre todolist donc!

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