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

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.

Mon premier test avec CartoDB

Dernièrement la firme Vizzuality de Madrid (maintenant aussi a New York) lance son programme CartoDB Beta pour donner un accès à PostGIS dans la Cloud via un API.  J'ai assisté à la présentation du produit au FOSS4G de Denver et je m'étais promis de l'essayer à sa sortie.   Justement hier soir suite aux élections en Espagne un des fondateurs de la firme (@saleiva) à fait ce petit démo avec cartodb pour visualiser les résultats 

Le produit est interessant et pourra aider très rapidement les développeurs désirant avec accès à des données géographiques via Postgresql / PostGIS.  En gros c'est un API genre Google fusion table  qui permet d'interagir avec une BD (select, insert, update, delete)  géo directement via un URL.  C'est assez bien fait à mon avis, et permettra d'éviter l'utilisation des services OGC plus lourd à implanter et supporter pour certaines catégories d'applications Web.

 

On a un problème de sécurité (SQL Injection) vous allez dire? Et oui c'est possible si vous ne protéger pas votre table de données.  Un billet d'un collègue Bills Dollins du Maryland a soulevé ce point le mois dernier et un développeur de CartoDB a répondu :

" ... the security model around CartoDB is based on public/private tables. When you make a table public what you are essentially doing is giving read permissions to a “public” user. That means that anybody can fire requests to your table, but they will hit PostgreSQL security if they try to modify something. So yeah, SQL injection as you want, as soon as you try to write, it will fail. You can actually try manually by taking the URL bill is calling on the example and run it directly."

C'est donc la responsabilité du programmeur d'inclure une sécurité via authentification "OAuth" pour protéger la table de données s'il le désire. Reste que l'API offre de belles options on peut par exemple visualiser rapidement les données d'une table via un petite app "embed_map"
On pourra aussi ajuster l'affichage des données on modifiant la thématique complexe via CartoCSS ou plus simplement via un contrôle manuel:
On aurra aussi quelques options pour configurer le fond de carte par défaut de votre API.  Pour l'instant on doit utiliser Google Maps mais j'imagine qu'on pourra éventuellement faire pointer la table sur le dépôt d'une carte tuilée:
Je n'ai pas fait une petite application pour tester l'API, mais j'ai fait quelques essais avec des requêtes.  Comme par exemple, un SQL "st_intersects" avec un Linestring en WKT.  Cette requête SQL retourne un segment dans le secteur du Vieux-Québec directement dans la carte configurée de mon compte:
Pour les amateurs de Geojson je dirais ceci:

Le potentiel est bon et on peut tester le service avec un maximum de 100 MB pour l'instant.  On peut imaginer qu'on aura un plan payable au mois, pour acheter plus d'espace ... On ne peut pas télécharger CartoDB et l'installer sur votre serveur.  À moins que le projet ne devienne Open Source mais rien ne transpire pour l'instant.  On peut télécharcher le projet sur Github (Merci a Alan Boudreault - mapgears.com).  En passant, le site pour charger et gérer les tables PostGIS est vraiment bien fait et est monté en Ruby on rails.

Charger un CSV avec psql dans Postgresql

Pour convertir un fichier CSV avec des coordonnées x,y en champ de type géométrie dans PostGIS, on peux utiliser la commande COPY en ligne de commande psql! 1) Copier le fichier CSV dans un répertoire visible du serveur Posrgresql. Dans le cas de Linux utiliser sudo et changer les droits du fichier CSV: sudo cp fichier_csv dans_repertoire_postgresql sudo chown postgres:postgres fichier_csv 2)Créer la table dans Postgres create table ma_table ( code character(10) NOT NULL, poste character(100) NOT NULL, adresse character(100) NOT NULL, municipalite character(100) NOT NULL, cp character(7) NOT NULL, _lat float, _long float, CONSTRAINT ma_table_pk PRIMARY KEY (code) ); 3)Charger les données du fichier CSV: set client_encoding='LATIN1'; COPY ma_table FROM '/srv/postgresql/scripts/ma_table.csv' delimiter ';'; 4)Ajouter la colonne geométrique select AddGeometryColumn ('public','ma_table','the_geom',4326,'POINT',2); 5)faire le update à partir des coordonnées x,y UPDATE ma_table set the_geom = ST_GeomFromText('POINT('||_long||' '|| _lat||')',4326); 6) Enjoy!

Un autre roadmap de création d’une BD spatiale sous Postgresql/PostGIS

Je reprends dans ce billet, la procédure complète pour monter rapidement une BD spatiale sous Postgresql/PostGIS. On utilisera souvent cette base de donnée simplement pour avoir un contenant de nos données géométriques.  Plusieurs utiliseront la chemin le plus court pour se connecter au serveur avec l’usager/mot de passe par défaut soit: postgres/postgres.  Il faut par contre idéalement, créer un autre utilisateur pour accéder au données SANS oublier de changer le mot de passe par défaut du serveur.

 

1) Installation du serveur Postgresql/PostGIS ET de l’utilitaire pgadmin pour un client ubuntu sudo apt-get update sudo apt-get install postgresql postgresql-client postgresql-contrib postgresql-8.4-postgis sudo apt-get install pgadmin3 sudo su postgres psql alter user postgres with password 'pw_defaut_a_changer_de_grace';

2) Un fois l’extension PostGIS intallée dans Postgresql, il ne reste qu’a trouver les deux fichiers SQL permettant d’accéder au fonctions spatiales et de la cherger avec psql.

Note: On va créer la BD template avec un encodage UTF-8.  Dépendant des versions à venir, le répertoire et même le nom des fichiers SQL peut varier... sudo su postgres createdb -E utf8 -O postgres postgis createlang plpgsql postgis psql -d postgis -f /usr/share/postgresql/8.4/contrib/postgis-1.5/postgis.sql psql -d postgis -f /usr/share/postgresql/8.4/contrib/postgis-1.5/spatial_ref_sys.sql 3) Vérifier l’installation Postgresql / PostGIS est bien correcte. On doit avoir la table geometry_table et spatial_ref_sys sudo su postgres psql postgis d Liste des relations Schéma |        Nom        | Type  | Propriétaire --------+-------------------+-------+-------------- public | geography_columns | vue   | postgres public | geometry_columns  | table | postgres public | spatial_ref_sys   | table | postgres 4) Par la suite on peux créer des BD (une bd par projet par exemple) avec le template postgis et un user sudo su postgres createdb -E utf8 -W -O postgres -T postgis osm psql -d osm CREATE USER osmusr WITH PASSWORD 'osm'; GRANT ALL PRIVILEGES ON DATABASE osm to osmusr; GRANT ALL ON TABLE geometry_columns TO osmusr; GRANT ALL ON TABLE spatial_ref_sys TO osmusr; alter user postgres with password 'petitpou'; q exit (sortir du super user postgres) psql -d osm -U osmusr -W