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.

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!

Scribe – Une nouvelle façon de faire du mapfile

Créer une belle carte qui véhicule efficacement son message n'est pas chose facile. On dispose souvent d'une grande quantité de données qu'on ne souhaite pas afficher en tout temps ou qu'on souhaite afficher mais avec un style qui varie en fonction du niveau d'échelle. Lorsqu'on défini le style d'une carte, chaque petit détail est important. Le processus se fait généralement de facon itérative et implique beaucoup de modifications, à tous les niveaux d'échelle.

Dans Mapserver, la gestion de l'échelle se fait avec les tags MINSCALEDENOM et MAXSCALEDENOM au niveau du LAYER ou de la CLASS. Cela veut dire qu'aussitôt qu'un élément du style change pour un LAYER ou une CLASS donné, à un certain niveau d'échelle, un nouveau LAYER ou une nouvelle CLASS doit être créé. On se trouve ainsi a réécrire beaucoup de texte pour ne modifier peut-être qu'un seul paramètre. La tâche devient encore plus fastidieuse si on doit changer un paramètre (par exemple la couleur d'une route) définie dans 16 LAYERs différents.

"Scribe" est un outil en python que nous avons créé pour faciliter l'écriture d'un "mapfile" en utilisant des variables et des raccourcis pour la gestion des échelles via des niveaux (maintenant la norme).  Cette façon de faire s'apparente à "Basemaps" mais est plus simple d'utilisation et généralement moins verbeuse.

Gestion des niveaux d'échelle

LAYER {
    1-16 {
        NAME: 'land'
        TYPE: POLYGON
        @layerconfig
        DATA {
            1-4: '110m_physical/ne_110m_land'
            5-10: '50m_physical/ne_50m_land'
            11-16: '10m_physical/ne_10m_land'
        }
        CLASS {
            STYLE {
                COLOR {
                    1-6: '#EEECDF'
                    7-16: '#AAA89B'
                }
                OUTLINECOLOR: 200 200 200
                OUTLINEWIDTH: @land_ol_width
            }
         }
     }
 }

Dans l'exemple précédent, un LAYER appelé "land" est créé. Le tag "1-16" signifie que ce LAYER doit etre affiché des niveaux d'échelle 1 à 16. Le script traduit automatiquement ces niveaux en MINSCALEDENOM et MAXSCALEDENOM. De plus, les données utilisées (DATA) des niveaux 1 à 4 sont au 110m tandis qu'elles sont au 50m des niveaux 5 à 10 et au 10m des niveaux 11 à 16. La couleur (COLOR) elle aussi change en fonction de l'échelle. En utilisant cette nomenclature, il devient très facile de modifier un paramètre pour un ou plusieurs niveaux d'échelle sans avoir à réécrire beaucoup de texte ou à faire des modifications à plusieurs endroits.

Définition et utilisation de variables

"Scribe" ne permet pas seulement de gérer les échelles mais il permet aussi de définir des variables réutilisables. Dans l'exemple ci-haut, les variables "layerconfig" et "land_ol_width" sont appelées a l'aide d'un "@". Ces variables sont définies de la facon suivante:

VARIABLES {
    layerconfig {
        GROUP: 'default'
        STATUS: ON
        PROJECTION {{
            'init=epsg:4326'
        }}
        PROCESSING: 'LABEL_NO_CLIP=ON'
        PROCESSING: 'CLOSE_CONNECTION=DEFER'
     }
     land_ol_width: 1
 }

La variable "layerconfig" contient plusieurs paramètres utilisés dans la définition de presque tous les LAYERs. Ainsi, pour chaque nouveau LAYER, il suffit d'écrire "@layerconfig" pour que tous les paramètres entrent dans la définition du LAYER. L'autre variable 'land_ol_width' prend une valeur unique.

Il est à noter que dans la définition de la variable "layerconfig", PROJECTION est suivi de deux "{". Cette syntaxe permet de gérer les tags comme PROJECTION, METADATA, PATTERN etc qui ne contiennent aucun paramètre, seulement du texte.

Blocs de commentaires

"Scribe" permet également d'utiliser des blocs de commentaire, ce qui n'est pas possible dans Mapserver seulement. D'autres types de commentaire sur une seule ligne peuvent aussi être écrits et apparaîtront dans le "mapfile" résultant.

LAYER {
    1-16 {
        NAME: 'land'
        TYPE: POLYGON
        @layerconfig
        DATA {
            1-4: '110m_physical/ne_110m_land'
            5-10: '50m_physical/ne_50m_land'
            11-16: '10m_physical/ne_10m_land'
        }
        CLASS {
            STYLE {
                COLOR {
                    1-6: '#EEECDF'
                    7-16: '#AAA89B'
                }                
                ##Les commentaires précédés par ## apparaissent
                ##dans le mapfile résultant.
                ##Les blocs de commentaires entre /* */
                ## n'apparaissent pas dans le mapfile résultant.
                /*
                OUTLINECOLOR: 200 200 200
                OUTLINEWIDTH: @land_ol_width
                */                
             }
         }
     }
 }

Conclusion

Finalement, pour exécuter le script, il suffit d'utiliser une seule commande, configurable avec certaines options:

python scribe.py

Le résultat est un mapfile parfaitement indenté avec gestion des échelles et souvent beaucoup plus de ligne que ce qu'il a fallu écrire. "Scribe" représente donc une économie importante en temps et produire des cartes de qualité devient plus facile et agréable.

Il est fort probable que "Scribe" évoluera et que de nouvelles fonctionnalités s'ajouteront.

Pour plus d'information ou pour utiliser "Scribe", consulter l'adresse suivante :

https://github.com/solutionsmapgears/Scribe.git Pour tout commentaire ou pour rapporter un bug : Charles-Éric Bourget cbourget@mapgears.com

L’échec d’Apple Maps

On avait bien hâte de la voir cette nouvelle Apple Maps.  Étant moi-même maniaque de cartographie Web, j'étais impatient de voir le résultat.  Les aperçus qu'on nous avait servis étaient prometteurs et j'espérais avoir un produit épuré, basé sur des données Open Street Map complété intelligemment avec d'autres sources de données.

L'appréciation générale de la nouvelle application de cartographie semble catastrophique. Selon Snappli (qui se spécialise dans la compression de données pour appareils mobiles), il ne resterait que 4% d'utilisateurs de Apple Maps. J'ai donc fait quelques petites validations dans la région de Québec.

La liste bien garnie de fournisseurs de données fut donc ma première surprise.  Règle générale au Québec, TomTom est le fournisseur officiel pour le réseau routier.  Pas vraiment surprenant, car pour la fonction de navigation routière ("turn-by-turn") il fallait obligatoirement utiliser TomTom.  Pour le reste on ne retrouve absolument rien de Open Street Maps pour les informations de contexte dans la carte!

Le principe fondamental "homogénéité" des données est pourtant élémentaire en cartographie.  L'utilisation de couches généralisées (hydrographie, boisée, etc) sans niveau de précision suffisant à cette échelle est une erreur de débutant, pour ne pas dire inacceptable!  La rivière Chaudière n'aboutit même pas jusqu'au fleuve, puisqu'on devrait gérer le réseau hydrique filamentaire à cette échelle.   Avec un peu de déception, je n'ai même pas poussé plus loin l'analyse, je vous laisse le soin de vous amuser les nombreux exemples d'erreurs répertoriées sur ce site qui deviennent pour le moins gênants pour le géant Apple.  Forcé d'admettre l'échec lamentable, c'est sans surprise que Tim Cook lui-même a du fournir des excuses aux utilisateurs iOS 6!

Apple a été tenté par des fonctionnalités sexy comme la navigation "turn-by-turn", le 3D, la reconnaissance vocale Siri, etc.   Mais on a oublié de se concentrer sur la base: une belle carte pratique à utiliser et pertinente!  La compagnie TomTom n'a pas été à la hauteur pour fournir ce dont Apple avait besoin.  Ils font de bons appareils de navigation automobile, mais TomTom n'a pas ce qu'il faut pour produire une cartographie Web de qualité pour les appareils mobiles.   De mon avis personnel, l'esthétique de la cartographie Web de TomTom dans la région de Québec est loin du minimum requis qu'une carte doit avoir de nos jours.  Je vous laisse le soin de voir vous-mêmes.

Investissements du Québec en santé

english version follow...
Nous avons réalisé dernièrement une nouvelle application de visualisation de données sur les "Investissements du Québec en santé".  Nous avons trouvé l'application vraiment "cool" alors on a décidé de la diffuser durant la présente campagne électorale 2012 parce qu'il est intéressant de pouvoir visualiser la taille de ces investissements.
Pour cette application, nous avons utilisé l'information publique du "Ministère de la Santé et des Services sociaux" diffusée via le site de données ouvertes du gouvernement du Québec, et des outils Web de visualisation de données pour mettre en relief le poids des investissements en santé versus la population et le budget total de la province de Québec.
Pour réaliser cette application, nous avons utilisé la librairie d3.js, Leaflet, puis Mapserver pour tuiler ( Mapcache ) et diffuser notre carte.
------------------------------
We recently launched our new data vizualisation application "Investissement du Québec en santé".  We find it cool so we decided to publish it.  It's interesting to visualize the size of these investments in this Quebec election campaign.
For this Web application, we used public data from the Government of Quebec Opendata Web site. We used data from "Ministery of Health and Social Services" and data visualization tools to highlight the amount of investments in health versus population of Quebec and the total budget of the province of Quebec.
To build this application, we used d3.jsLeaflet, and we used Mapserver to tiled (Mapcache) and served our map.

Solutions Mapgears s’installe à Québec

Après un début d'année passablement actif pour Solutions Mapgears, le temps était venu de s'installer dans la région.   C'est donc un grand plaisir que nous pouvons annoncer notre nouvelle place d'affaires dans la ville de Québec au 2383 chemin Sainte-Foy, bureau 202.

De plus, comme une bonne nouvelle ne vient jamais seule,  c'est avec une grande fierté que nous avons récemment accueilli dans l'équipe du groupe Mapgears, deux finissants en génie géomatique de l'université Laval Charles-Éric Bourget et Marc-André Barbeau.

Le groupe Mapgears poursuit sa croissance et garde le cap sur ses objectifs à long terme.  L'équipe compte maintenant 9 employés (3 à Québec et 6 à Chicoutimi) et continue d'aider ses clients tant à l'échelle locale, nationale et aussi internationale qui représentent maintenant près de 50% de son chiffre d'affaire.

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.

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/  

Un WMS de AQréseau

J'ai fait quelques essais avec la couche de données gratuite AQréseau. J'ai chargé les données dans Postgresql/PostGIS. Ce WMS (http://mapcoop.org/cgi-bin/aqreseau) diffusé avec Mapserver est disponible pour consultation dans votre SIG bureautique préféré. La cartographie du service n'est pas parfaite mais elle montre ce qu'il est possible de faire avec cette nouvelle couche de données. Je l'ai aussi placé dans une vue OpenLayers pour rendre la chose plus facile et la consulter rapidement. Le WMS AQréseau: http://mapcoop.org/cgi-bin/aqreseau Du point de vue purement cartographique, on pourra y faire quelques découvertes.  Entre autres la liste de données qu'il est possible de cartographier. Le nombre de segments "Pistes cyclables" (11) n'est pas très élevé pour l'instant. Même chose pour les "Sentiers piétonniers"(2). Espérons qu'on aura accès à cette information prochainement... J'ajoute quelques statistiques sur la classification des données pour cette première version. select clsrte,caractrte,count(gid) from aq_routes group by clsrte,caractrte order by clsrte,caractrte;
clsrtecaractrteNb de segments
Accès localités isolées178
Accès ressourcesBretelle3
Accès ressources694
ArtèreBretelle263
ArtèreTraverse530
Artère19485
AutorouteBretelle4772
AutorouteTraverse301
AutorouteVirage en U1325
AutorouteVoie de service528
Autoroute8179
Collectrice564
Collectrice de transitBretelle63
Collectrice de transitCarrefour giratoire6
Collectrice de transitTraverse28
Collectrice de transitVirage en U14
Collectrice de transit15176
Collectrice municipaleBretelle309
Collectrice municipaleCarrefour giratoire14
Collectrice municipaleConflit SDA4
Collectrice municipaleTraverse799
Collectrice municipaleVirage en U5
Collectrice municipale34485
LocaleBretelle585
LocaleCarrefour giratoire12
LocaleConflit SDA30
LocaleDivers 51
LocaleLiaison maritime42
LocaleTournebride1
LocaleTraverse3168
LocaleVirage en U55
LocaleVirtuel617
LocaleVoie de service4
Locale280914
NationaleBretelle675
NationaleCarrefour giratoire78
NationaleTournebride15
NationaleTraverse131
NationaleVirage en U157
NationaleVoie de service9
Nationale20170
Piste cyclable11
RégionaleBretelle172
RégionaleCarrefour giratoire29
RégionaleTournebride5
RégionaleTraverse70
RégionaleVirage en U23
Régionale13561
Rue piétonne14
Sentier piétonnier2