#Pgrouting: Analyse comparative d’Isodistances de 500m à partir de Natntes métropole.

isochrones.jpgEn cartographie, une courbe isochrone est une courbe géométrique délimitant les points accessibles par un véhicule – terrestre ou aérien – en un temps donné (par exemple, la zone pouvant être desservie en moins de 30 minutes par un livreur de pizza ou un dépanneur de matériel informatique).Ou bien en distance parcourue dans ce cas précis on parle d’Isodistance.

1.      Objectif

Analyse comparative d’Isodistances de 500m à depuis de Natntes métropole.

2.      Méthodes de création d’isochrones

D’abord, pour créer une structure topologique complète et vous familiariser avec pgrouting, je vous  renvoie  à l’article précédente dans cette page: https://diouck.wordpress.com/postgispgrouting/

Après avoir modéliser notre réseau, nous allons calculer et visualiser les nœuds et arcs qui se situent à 500m de Nantes métropole. Pour ce faire nous allons utiliser l’algorithme driving distance de pgRouting pour calculer la distance la plus coute accessible à moins de 500m.

ØDriving distance pour les nœuds
--Drop table routing.pgr_drivingdistance_pt;
create table routing_nantes.pgr_drivingdistance_pt as
SELECT id1 AS node_id , cost,geom
FROM pgr_drivingdistance(
'SELECT edge_id as id, start_node as source, end_node as target, tps_distance as cost from routing.edge_data',
12616, 0.5, false, false
) as di
JOIN routing.node pt
ON di.id1 = pt.node_id;

Ø Driving distance pour les tronçons
--Drop table routing.pgr_drivingdistance_lgn;
create table routing_nantes.pgr_drivingdistance_lgn as
SELECT id1 AS node_id ,pt.tps_distance as cost,geom
FROM pgr_drivingdistance(
'SELECT edge_id as id, start_node as source, end_node as target, tps_distance as cost from routing.edge_data',
12616, 0.5, false, false
) as di
JOIN routing.edge_data pt
ON di.id1 = pt.start_node;   

arcs et noeuds avec pgr_drivingdistance.jpg

Comme vous pouvez le constater les deux résultats sont légèrement différents.

Les nœuds retournent exactement la distance entre 0 et 499m
tandis que les arcs(tronçons de voirie) dépassent de plusieurs mètres sur les extrémités.

Cette erreur est liée aux intersections entre les nœuds extrêmes et les arcs qui les touches au-delà des 500m. Maintenant nous allons comparer différents isodistances  à partir des nœuds et des arcs et voir lesquels s’approchent le mieux à la réalité.

a- Isodistance avec la fonction pgr_alphAShape

Cette fonction retourne un tableau avec des lignes (x, y) qui décrivent les sommets d’une forme alpha. Nous allons donc le coupler avec l’algorithme de driving distance pour créer un polygone à partir des sommets des nœuds et des arcs. Pour plus d’information je vous renvoie vers cette documentation:
A Closer Look at Alpha Shapes in pgRouting

--DROP TABLE routing.isochrone_pgr_alphAShape
CREATE TABLE routing.isochrone_pgr_alphAShape AS
SELECT ST_MakePolygon(ST_AddPoint(foo.openline, ST_StartPoint(foo.openline)))::geometry, 2 as alphaPoly
from (select st_makeline(points order by id) as openline from
(SELECT st_makepoint(x,y) as points ,row_number() over() AS id
FROM pgr_alphAShape('SELECT node_id::integer as id, st_x(geom)::float as x, st_y(geom)::float as y FROM routing.pgr_drivingdistance_pt')
) as a) as foo;

isochrone_pgr_alphAShape.jpg

L’Isodistance regroupe bien les sommets des nœuds pour créer un polygone. Cependant on peut constater une légère surestimation des distances.

b- Isodistance avec la fonction pgr_pointsAsPolygon

Comme la fonction précédente, elle va créer un polygone plus précis au tour des nœuds.

--DROP TABLE routing.pgr_pointsAsPolygon;
CREATE TABLE routing.pgr_pointsAsPolygon AS
SELECT 500, ST_SetSRID(geom,2154)
FROM
pgr_pointsAsPolygon(
'SELECT node_id::integer as id, st_x(geom)::float as x, st_y(geom)::float as y FROM routing.pgr_drivingdistance_pt') as geom;

On Remarque une très grande précision sur les noeuds. Mais il sous-estime l’emprise des noeuds. Essayons maintenant de fermer le polygone ouvert à l’intérieur.

--DROP TABLE routing.pgr_pointsAsPolygon_ST_ExteriorRing;
CREATE TABLE routing.pgr_pointsAsPolygon_ST_ExteriorRing AS
with pgr_pointsAsPolygon as (
SELECT 500, ST_SetSRID(geom,2154) as geom
FROM
pgr_pointsAsPolygon(
'SELECT node_id::integer as id, st_x(geom)::float as x, st_y(geom)::float as y FROM routing.pgr_drivingdistance_pt') as geom
)
SELECT ST_MakePolygon(ST_ExteriorRing(ST_GeometryN(geom,1))) as geom FROM pgr_pointsAsPolygon;

pgr_pointsAsPolygon_ST_ExteriorRing.jpg

 

c- Isochrone avec St_buffer sur les coûts (distance)

L’objectif est de pondérer les coûts en fonction de la distance. Cette méthode est inspirée de cet article http://anitagraser.com/2013/07/07/public-transport-isochrones-with-pgrouting/

--Isochrones à partir de buffer
drop table routing.buffer_cost;
CREATE TABLE routing.buffer_cost as
with buffer as (
select
case when cost<1 then (0.5-cost)*200 end as cost,geom
from routing.pgr_drivingdistance_pt)
select st_union(st_buffer(geom,cost)) as geom from buffer ;

isochrone_st_buffer_cost.jpgIl faudra adapter la requête. La méthode est bonne mais on remarque qu’en même des discontinuités spatiales. Nous avons des zones discontinues et inaccessibles alors qu’en réalité elles devraient l’être.

d- Isochrone avec st_concavehull et st_union

Cette méthode va nous permettre de créer une enveloppe concave autour d’un nuage de points (concave hull). Pour plus d’information je vous renvoie vers la documentation compète :
http://www.portailsig.org/content/sur-la-creation-des-enveloppes-concaves-concave-hull-et-les-divers-moyens-d-y-parvenir-forme

drop table routing. st_concavehull_st_union;
CREATE TABLE routing. st_concavehull_st_union as
SELECT 1 AS id, st_concavehull(st_union(t.geom), 0.7)--ST_ConcaveHull(ST_Collect(geom),0.95,true) --
FROM
(SELECT seq, id1 AS node, id2 AS edge, cost, pt.geom
FROM pgr_drivingdistance(
'SELECT edge_id as id, start_node as source, end_node as target, tps_distance as cost from routing.edge_data',
12616, 0.5, false, false
) AS di
JOIN routing.node pt
ON di.id1 = pt.node_id) t;

isochrone_st_concavehull_st_union.jpg

Vous pouvez modifier la précision à votre guise. Ici elle est de 0.7. L’isochrone reste correct mais imprécis avec quelques imperfections.

e- Isochrone avec St_ ConcaveHull et st_Collect

La méthode inspirée de ce poste :
http://gis.stackexchange.com/questions/95771/creating-isochrones-in-postgis-osm2po-pgrouting-and-then-saving-the-isochrones-a
La tolérance optimale se situe à 0.7. Vous pouvez évidemment l’adapter à vos besoins.

 Avec les arcs
drop table routing.ST_ConcaveHull_ST_Collect_lgn;
create table routing.ST_ConcaveHull_ST_Collect_lgn as
SELECT ST_ConcaveHull(ST_Collect(geom),0.70,false)
FROM routing.edge_data
JOIN (SELECT * FROM pgr_drivingdistance('
SELECT edge_id as id, start_node as source, end_node as target, tps_distance as cost from routing.edge_data',
12616, 0.5, false, false ))
AS route
ON routing.edge_data.start_node = route.id1 ;
 Avec les noeuds
drop table routing.ST_ConcaveHull_ST_Collect_pt;
create table routing.ST_ConcaveHull_ST_Collect_pt as
SELECT ST_ConcaveHull(ST_Collect(geom),0.70,false)
FROM routing.edge_data
JOIN (SELECT * FROM pgr_drivingdistance('
SELECT edge_id as id, start_node as source, end_node as target, tps_distance as cost from routing.edge_data',
12616, 0.5, false, false ))
AS route
ON routing.edge_data.start_node = route.id1 ;

Isochrones avec St_ConcaveHull et st_Collect.jpg

On a soit une surestimation liée aux intersections avec les nœuds externes sur les arcs de plus 500m ou une sous-estimation du polygone créé.

 

f- Isochrone avec ST_MakePolygon et ST_ExteriorRing

Objectif: nn ferme les polygones après un buffer pour créer un polygone unique.

 Avec les arcs
drop table routing.ST_MakePolygon_ST_ExteriorRing;
create table routing.ST_MakePolygon_ST_ExteriorRing as
WITH buffer_itineraire as (
SELECT --id1 AS node_id ,pt.cost,
st_buffer(st_union(geom),10 ) as geom
FROM pgr_drivingdistance(
'SELECT edge_id as id, start_node as source, end_node as target, tps_distance as cost from routing.edge_data',
12616, 0.5, false, false
) as di
JOIN routing.edge_data pt
ON di.id1 = pt.start_node)
---Fermer les polygones
SELECT ST_MakePolygon(ST_ExteriorRing(ST_GeometryN(geom,1))) as geom FROM buffer_itineraire;


 Avec les arcs et les noeuds
drop table if exists routing. ST_MakePolygon_ST_ExteriorRing_v2;
create table routing. ST_MakePolygon_ST_ExteriorRing_v2 as
with buffer_itineraire as (
SELECT et.edge_id ,
st_buffer(et. geom,10, 'endcap=square join=round') as geom
, 1 as factor
FROM
(SELECT id1,cost from pgr_drivingDistance(
'SELECT edge_id as id,start_node as source , end_node as target, tps_distance as cost FROM routing.edge_data',
'12616','0.5',false,false)
) firstPath
CROSS JOIN
(SELECT id1,cost from pgr_drivingDistance(
'SELECT edge_id as id,start_node as source , end_node as target,tps_distance as cost FROM routing.edge_data',
'12616','0.5',false,false)
) secondPath
INNER JOIN routing.edge_data et
ON firstPath.id1 = et.start_node
AND secondPath.id1 = et.end_node)
SELECT ST_MakePolygon(ST_ExteriorRing(ST_GeometryN(st_union(geom),1))) as geom FROM buffer_itineraire;

Isochrone avec ST_MakePolygon et ST_ExteriorRing.jpg

Ce sont les polygones qui s’approchent mieux de la réalité et en particulier l’isochrone avec les nœuds.

 

g- Synthèse

L’objectif était de voir comment créer un isochrones. Il existe comme vous pouvez le constater plusieurs façons d’y arriver. Les méthodes sont à adapter en fonction de vos besoins. L’isochrone qui se rapproche le plus de la réalité dans cette exercice est celle que j’ai optimisé avec un ST_MakePolygon couplée de ST_ExteriorRing sur les nœuds et non sur les troncons du réseau de voirie de Nantes métropole.

A suivre prochaine étapes les isochrones sachant qu’on peut opérer de la même façon pour comparer l’accessibilité d’un lieu  en minutes ou heures.

isochrones.jpg

 

 

 

#Monalisa2PostGIS: Importation de nombreux raster (Traitement par lots)

monalisa

Pour ceux qui travaillent avec Postgis de manière plus ou moins fréquente ,  si vous avez des difficultés d’intégration de plusieurs raster en même temps, voici un  script batch  qui fera office d’ETL. Évidemment les applications sont multiples par exemple regrouper plusieurs dalles d’ortho ou d’images satellites et calculer des indices de végétation.

Le scripte peut être exécuté sans problème malgré un nombre important de raster ( Tester sur un dossier de 50 dalles de 1km).  Voici les paramètres à modifier:

  1- D’abord définir les paramètres de connections

REM Parametre de connections
set PGPORT=5432 
set PGHOST=localhost
set PGUSER=postgres
set PGPASSWORD=postgres

2- Définir les  exécutables de postgres

REM  chercher le programme 
cd "C:\Program Files\PostgreSQL\9.5\bin"

3- Transformer en SQL l’image dans le schéma raster

Monalisa est disponible sur sa page  wikipedia

raster2pgsql   -d -I -C -e -Y   -s 4326   -t 128x128    "C:\raster\monalisa\*.jpg"  raster.monalisa >    "C:\raster\monalisa.sql" 
Avec -d pour supprimer la table si elle existe
-C pour la créer si elle n'existe pas
-I pour créer un overview
-s pour la projection et -t pour le tuilage

4- Importer le SQL dans la base de donnee bdu

psql -d bdu -f "C:\raster\monalisa.sql"

Pour info le script bat est disponible ici

5- Et visualiser votre résultat sur Qgis.

 

 

 

 

 

Mise en place d’un système de routage et de calcul d’itinéraires intégrable à PostgreSQL/PostGIS : « Pgrouting »

Dijkstra voiture direction 2.png

 

Il existe de nombreuses recherches sur l’utilisation des algorithmes de pgRouting pour l’analyse du plus court chemin et le calcul des isochrones. PgRouting est une extension de PostgreSQL et PostGIS. Il fournit divers algorithmes comme All Pairs Short Path, Bi-directional Dijkstra, and A*, Shortest Path y compris Turn Restriction, Driving Distance etc… et d’autres fonctions d’analyses de graphes.

Les fonctions pgRouting sont basées sur les coûts, qui sont dynamiques par opposition coûts de pré-calculés et peuvent donc s’appliquer à différents types de scénarios réels sur un réseau routier.

En dépit de la diversité des possibilités offertes par pgRouting pour l’analyse d’un réseau, et malgré d’intenses recherches sur la toile, il existe peu d’exemples sur la réalisation d’isochrones parfaites ou sur la gestion multimodale des coûts (vitesses/temps) sur un réseau.

La plupart des tutoriels ou des recherches existantes se concentrent essentiellement sur sa facilité d’utilisation, sur la recherche du plus court chemin qui elle-même est basé sur la notion de coût.

Dans ce manuel, nous allons apprendre les différentes étapes de modélisation d’un réseau routier. Pour ce faire j’ai fait le choix de travailler uniquement avec des données gratuites d’OpenStreetMap au lieu des données d’IGN, TéléAtlas et autres qui sont payantes.

Nous allons donc voir aussi comment installer pgRouting, importer les données d’OpenStreetMap , créer un graphe planaire orienté et l’exploiter dans Qgis. Puis nous allons modéliser les coûts (voiture, peton, vélo et distance) en respectant le sens de la circulation. Je vous montrerai également comment créer des zones de chalandises ou isochrones.

Tutoriel accessible ici:

 

 

#pgRouting : Schéma Directeur d’Accessibilité voirie et espace public de Nantes Métropole (SDA)

1.      Cadre de la demande

La politique menée par Nantes Métropole en faveur de l’accessibilité depuis de nombreuses années.  La volonté d’aller au-delà du cadre réglementaire dans la construction d’un territoire accessible pour toutes et tous) en favorisant au maximum la continuité de la chaine de déplacement.  La réalisation d’un document cadre stratégique issu de l’expertise urbaine et de la concertation des usagers, permettant d’objectiver les décisions. La constitution, à terme, d’un SIG accessibilité, guide pour les futurs travaux d’aménagements de la Métropole.

2.      Objectifs

Identifier les itinéraires susceptibles d’être utilisés de façon prioritaire par des personnes à mobilité réduite (et de manière générale par les piétons) à travers la ville, du domicile au lieu de travail, d’achat, de loisirs, …Définir le volume des cheminements à diagnostiquer sur le Territoire de Nantes Métropole (analyse du niveau d’accessibilité des itinéraires identifiés comme prioritaires)Aider à la priorisation de la mise en accessibilité des arrêts TC prévue dans le cadre de la Stratégie du SDAp’TC 2015-2020.

3.      Méthodologie

Application avec pgRouting

Cheminements sur voirie, de moins de 300 mètres, reliant les Etablissements Recevant du Public (ERP) aux stationnements réservés aux personnes à mobilité réduite ou aux arrêts de transport en commun les plus proches, ainsi que l’intégralité des Plans d’Accessibilité de la Voirie et de l’Espace Public (PAVE) élaborés par les communes

Schéma directeur d'Accessibilité de Nantes Métropole (SDA).png

Pour en savoir plus l’étude est accessible ici

Rapport dynamique: Liste des attributs d’une table dans #Postgre

Disposant d’une base de données GeoNework , je me suis interrogé sur comment intégrer la liste des attributs  d’une table dans  sa fiche de métadonnées. Et surtout comment les visualiser d’une manière dynamique.  Alors j’ai regardé du coté  PGAdmin, et oui on peut générer directement un rapport en html. Mais comment le rendre dynamique?

D’abord après de fructueuses recherches voici la requête qui me permet d’afficher pour une table donnée:  les attributs, les types, le nom de la table, les commentaires, s’il est rempli ou pas par défaut, la clef primaire…

Voici un exemple avec la table bâtiment dans le schéma majic (données foncières)


 SELECT a.attnum as numero,
 a.attrelid::regclass as table
 ,a.attname AS attribut
 ,format_type(a.atttypid, a.atttypmod) AS type
 ,a.attnotnull AS notnull
 ,coalesce(p.indisprimary, FALSE) AS primary_key
 ,f.adsrc AS default_val
 ,d.description AS commenatire
 FROM pg_attribute a
 LEFT JOIN pg_index p ON p.indrelid = a.attrelid AND a.attnum = ANY(p.indkey)
 LEFT JOIN pg_description d ON d.objoid = a.attrelid AND d.objsubid = a.attnum
 LEFT JOIN pg_attrdef f ON f.adrelid = a.attrelid AND f.adnum = a.attnum
 WHERE a.attnum > 0
 AND NOT a.attisdropped
 AND a.attrelid = 'majic_2014.d44_2014_batiment'::regclass -- A adapter
 ORDER BY a.attnum;
 

Voici le résultat:

majic.png

Et avec le rapport de pgAdmin(fichier/rapport rapide) ça donne ça:

rapport pgadmin.png

 

Alors tout ça c’est magique. Mais comment le rendre dynamique pour éviter de générer des milliers de pages tml? Je suis donc passé par du php pour transmettre des données par url.

Par exemple,  vous voulez afficher les attributs de la table « majic_2014.d44_2014_lotslocaux »   voici l’url à transmettre dans GeoNework :

http://localhost/developpement/rapport/rapport.php?table=majic_2014.d44_2014_lotslocaux
Vous remarquerez que la variable "table" permet de récupérer les attributs de la table que vous voulez afficher

Voici le résultat dans GeoNework :metadonnees.png

Voici le résultat :

rapport postgresql.png

Vous trouverez le code ici :

Évidemment vous pouvez également l’utiliser pour générer directement d’autres requêtes.

Pour plus de détails concernant la transmissions de variables ici

 

 

 

Calcul de l’indice de végétation NDVI avec PostgeSQL/Postgis

En télédétection, les indices font parties des méthodes de traitement que l’on appelle les transformations multispectrales. Ils consistent à convertir les luminances mesurées au niveau du capteur satellitaire en grandeurs ayant une signification dans le domaine de l’environnement.

Basés sur le caractère multispectral des données satellitaires, ils permettent de décrire l’état d’un phénomène. Un indice de végétation par exemple, peut rendre compte du stade de croissance végétale à un moment donné.
Tous les indices, que ce soient les indices de végétation, les indices des sols, les indices relatifs à la colonne d’eau, etc., reposent sur une approche empirique basée sur des données expérimentales. Les indices de végétation sont très utilisés d’une part, pour identifier et suivre la dynamique de la végétation, mais aussi pour estimer certains paramètres biophysiques caractéristiques des couverts végétaux, comme la biomasse, l’indice de surface foliaire, la fraction de rayonnement photosynthétique actif, etc.

L’indice le plus connu et le plus utilisé est l’indice de végétation par différence normalisé ou indice de Tucker (NDVI en anglais) (Rouse and Haas, 1973 ; Tucker, 1979). Son expression est la suivante :

(pir-rouge)/(pir+rouge)

Je suppose que vous vous êtes déjà familiarisés avec Postgis ,Qgis et Gdal.

Les données : Images Quickbird (Visible +proche infrarouge)

composition colorée DIOUCK

Composition colorée Rouge +PIR

Assemblage des données avec Gdal

Pour faciliter les traitements, j’ai assemblé les bandes du rouge et du proche infra-rouge dont nous aurons besoin pour calculer le NDVI. Si vous n’avez pasGdal je vous conseille d’intaller FWtools

–Stack Bande
gdal_merge -separate  bande3.tif bande4.tif -o bande34.tif

Pour vérifier notre stack voici la commande gdalinfo bande34.tif

gdalinfo_bande

 gdalinfo sur le stack

Et là on distingue bien nos deux bandes : La bande 1 correspond au rouge et la bande 2 au PIR

Importation des données dans Postgis :

J’ai créé un script bat pour l’intégration des données

REM Paramètres de connections
set PGPORT=5432
set PGHOST=localhost
set PGUSER=postgres
set PGPASSWORD=postgres
REM on traite le raster en le transformant en une suite de commandes SQL pour le schema ‘raster’ . Le –t permet de tuiler la donnée et le –I une indexation pour  gagner en performance
« C:\Program Files\PostgreSQL\9.3\bin\raster2pgsql »  -s 4326 -C -I -r -M -F -t 100×100  « D:\Postgis Raster\Satellite\QuickBird\ bande34.tif »  raster. bande34> bande34.sql
REM on crée la table raster dans notre base BDU de donnée via psql
« C:\Program Files\PostgreSQL\9.3\bin\psql » -d bdu -f bande34.sql
 
 
 
 

Calcul de l’indice de végétation  NDVI

L’indice le plus connu et le plus utilisé est l’indice de végétation par différence normalisé ou indice de Tucker (NDVI en anglais) (Rouse and Haas, 1973 ; Tucker, 1979). Son expression est la suivante :

La formule est la suivante (pir-rouge)/(pir+rouge)

Voici la requête dans postgis

create table ndvi as
select ST_MapAlgebra(a.rast, 2, a.rast, 1, ‘([rast1] – [rast2])/([rast1] + [rast2])::float’, ’32BF’) as rast
FROM raster.bande34 a

 

Visualisation dans QGIS

Comme vous pouvez le constater nous avons bien les valeurs du NDVI comprises entre -0.9 à 0.68.
Une simple classification permet de visualiser les différentes classes d’occupation du sol qui sont pareilles que la classification que j’ai fait avec les mêmes données à partir de QGIS

NDVI avec POSTGIS

NDVI POSTGIS DIOUCK
NDVI avec QGIS
NDVI QGIS DIOUCK

 

 

Système de référence linéaire

Le système de référence linéaire (LRS) est un système de référence qui permet de localiser les éléments à l’aide d’une mesure le long d’un élément linéaire. Chaque entité est localisée par un point connu sous le nom de « jalon » ou par un événement linéaire (« segment »). Le système est conçu de telle sorte que si un segment de route est modifié, seuls les jalons du segment modifié doivent être mis à jour.

Un LRS est approprié pour la gestion des données relatives aux entités linéaires comme les routes, les chemins de fer, les pipelines (de transport de pétrole et gaz), et des rivières.
Un LRS permet d’identifier l’emplacement des entités et des caractéristiques d’un pipeline par mesure de distance depuis le début de l’oléoduc. Un LRS est fourni par exemple par Arcgis, GRASS et PostGIS.

Objectif :
Dans le cadre de la mise en place d’une application web mapping destinée aux travaux et à la circulation dans Rennes Métropole, le service SIG en étroite collaboration avec la DIRO va produite un référentiel linéaire pour cartographier le trafic en temps réel. Après une étude approfondie de la bibliographie nous nous sommes orientés sur le référencement linéaire de Postgis.

Outils méthodologiques :

Les données :

Système de référence linéaire

Système de référence linéaire

Nous disposons les données référentielles de la DIRO :
Liste des tronçons : Un fichier qui nous donne une information permettant de tracer les tronçons à partir des points de références.( troncons_tab.csv)
Points de référence : Cette donnée contient 2773 points de références sur l’ensemble de la Bretagne.
(BORNAGE_ROUTES_DIRO_2012.shp)
Et les routes : Routes départementales et nationales sur l’ensemble de la Bretagne (route_diro_2012shp).

Fonctionnalité du référencement linéaire de Postgis :
ST_line_locate_point(LineString, Point): Opérateur qui permet de localiser en pourcentage un point sur une ligne avec en entrée la géométrie de la ligne et celle du point.
ST_line_interpolate_point(linestring, location) :Cet opérateur permet de créer projeter un point sur une ligne en créant un point
ST_line_substring(linestring, start, end): Permet de tracer une ligne à partir des points(début et fin) sur une ligne quelconque.
Exemple d’application : Tronçon= « 8005 » sur la N136
Nous allons localiser les points de références avant de les reprojeter sur les routes. Et enfin tracer le troncons.

-- DROP TABLE a84_8006lgn;
CREATE TABLE route_8005_lng as
SELECT tr.id,tr.route,sens,
st_line_substring(CASE
WHEN geometrytype(rt.the_geom) = 'MULTILINESTRING'::text THEN geometryn(rt.the_geom, 1)
ELSE rt.the_geom
END,m_pr_debut,m_pr_fin),
rt.gid + 1 AS objectid
FROM troncons_tab tr, "BORNAGE_ROUTES_DIRO_2012" bor, route_diro_2012 rt,

(select (((length(rt.the_geom)+ split_part(tr.pr_debut::text, ‘+’::text, 2)::double precision) * st_line_locate_point(
CASE
WHEN geometrytype(rt.the_geom) = ‘MULTILINESTRING’::text THEN geometryn(rt.the_geom, 1)
ELSE rt.the_geom
END, bor.the_geom) + split_part(tr.pr_debut::text, ‘+’::text, 2)::double precision) / (length(rt.the_geom)+ split_part(tr.pr_debut::text, ‘+’::text, 2)::double precision)) as m_pr_debut
from troncons_tab tr, pr_intersect_route bor, route_diro_2012 rt
WHERE tr.id::text = ‘8005’::text AND split_part(tr.pr_debut::text, ‘+’::text, 1) = bor.pr::text AND st_intersects(rt.the_geom, bor.the_geom) AND rt.gid = 15
) AS prdeb,

(select (((length(rt.the_geom)+ split_part(tr.pr_fin::text, ‘+’::text, 2)::double precision) * st_line_locate_point(
CASE
WHEN geometrytype(rt.the_geom) = ‘MULTILINESTRING’::text THEN geometryn(rt.the_geom, 1)
ELSE rt.the_geom
END, bor.the_geom) + split_part(tr.pr_fin::text, ‘+’::text, 2)::double precision) / (length(rt.the_geom)+ split_part(tr.pr_fin::text, ‘+’::text, 2)::double precision)) as m_pr_fin, tr.id
from troncons_tab tr,pr_intersect_route bor, route_diro_2012 rt
WHERE tr.id::text = ‘8005’::text AND split_part(tr.pr_fin::text, ‘+’::text, 1) = bor.pr::text AND st_intersects(rt.the_geom, bor.the_geom) AND rt.gid = 15

) AS prfin,

( select CASE
WHEN split_part(tr.pr_fin::text, ‘+’::text, 2)::double precision < split_part(tr.pr_fin::text, ‘+’::text, 2)::double precision THEN ‘PRcroissant’
ELSE ‘PRdecroissant’
END as sens
from troncons_tab tr
where tr.id::text = ‘8005’::text) as sens
WHERE tr.id::text = ‘8005’::text AND split_part(tr.pr_debut::text, ‘+’::text, 1) = bor. »PR »::text AND st_intersects(rt.the_geom, bor.the_geom) AND rt.gid = 15;
ALTER TABLE route_8005_lng OWNER TO geotravaux;

Resultat

Système de référence linéaire

Système de référence linéaire