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.

2 comments — post a comment

smercier

Cool Martin, oui jva faire ça, sauf les commentaires sont en français :-( tk jva poster sur ton github et peut-être le traduire si ya du temps
cheers

Leave a Reply

Your email address will not be published. Required fields are marked *

*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

Current month ye@r day *