#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

 

 

 

Accessibilité des équipements & loisirs à Paris

Un site informatif responsive design conçu avec les données de l’Open data de la ville de Paris. Vous y trouverez uniquement les équipements gérés par la Ville de Paris pour lesquels la Ville a la responsabilité de la mise en accessibilité, par exemple uniquement les crèches municipales et non les crèches associatives, ou uniquement les lycées professionnels, gérés par la Ville, et non les lycées d’enseignement général, gérés par la Région.

Dans un premier temps, seuls sont traités dans cette carte les handicaps moteurs, auditifs et visuels, le handicap mental fera l’objet d’un recensement ultérieur. De même les services travaillent à l’accessibilité de la carte elle-même.

Sont définis comme accessibles les établissements où vous pouvez bénéficier du service proposé par l’établissement ou l’équipement, même si l’accès à la totalité des locaux n’est pas toujours possible.

Niveaux d’accessibilités

Non accessible
L accueil du public est accessible pour tous
Tous les services au public sont accessibles pour tous
Tous les locaux ouverts au public sont accessibles pour tous
Niveau 3 ERP + tous les locaux de travail sont accessibles pour le personnel
Donnée non renseignée

Accessibilité des équipement de Paris »:
http://www.diouck.esy.es/accessibilite/accessibilite.html

accessibilite

accessibilité

 

Responsive design mobile

Responsive design mobile

#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.

 

 

 

 

 

#Postgres / #PostGIS: Importation de nombreux shapefiles (Traitement par lots)

       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 shapefiles en même temps, voici un  script batch  qui fera office d’ETL.

     Il peut être exécuté sans problème malgré un nombre important de shapefiles ( Tester sur un dossier de 55 shp ).  Voici les paramètres à modifier:

  • D’abord définir les variables vers les shp, et exécutables de postgres
set PATH_SHP="C:\Users\ABDOU\shapefiles"
set PATH_PSQL=C:\Program Files\PostgreSQL\9.5\bin\psql.exe
set PATH_SHP2SQL=C:\Program Files\PostgreSQL\9.5\bin\shp2pgsql.exe
set PATH_SQL="C:\Script_Abdou\shp\shp2pgsql\sql"
  • Les variables pour les paramètres de connexions
set pghost=votre host exemple localhost
set pgport=5432
set pgdb=votre base de données
set pgsrid= votre système de projection exemple 2154
set pggeom= le nom de votre colonne géométrique exemple geom
set pgencoding="CP1252"
set pgschema= le schema ou importer les shp
  • La procédure d’exécution pour chaque shp
REM "Parcourir le dossier et créer un fichier sql pour chaque shp "
for %%f in (%PATH_SHP%\*.shp) do "%PATH_SHP2SQL%" -s %pgsrid%   -g %pggeom% -d -D -i -I -W %pgencoding% %%f %pgschema%.%%~nf > %PATH_SQL%\%%~nf.sql

REM "Parcourir le dossier sql et intégrer les fichiers SQL dans la base de données Postgres et créer un fichier log pour chaque données intégrées"
for %%f in (%PATH_SQL%\*.sql) do "%PATH_PSQL%" -h %pghost% -p %pgport% -d %pgdb% -L %%~nxf.log -U %pguser% -f %%f

REM " Supprimer tous les fichiers sql dans le dossier SQL "
for %%f in (%PATH_SQL%\*.sql) do del %%f

Après ces changements enregistrer le fichier au format .bat
Exemple  shapefiles2postgis.bat
Double cliquer sur le fichier .bat pour l’exécuter
Et attendre les résultats ….🙂

Pour plus de détails le fichier est disponible ici

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:

 

 

Le marché de l’immobilier de bureaux à Nantes Métropole, de 2007 à 2015

onita2015

Cette carte présente l’évolution des surfaces transactées sur le marché de l’immobilier tertiaire à Nantes Métropole, depuis 2007.
Chaque opération a été localisée sur la carte, et présente la surface transactée (en m2), l’état des bureaux (neuf ou seconde main) et la nature de la transaction (vente ou location).

Application accessible ici

Appli web développeé avec leaflet et inspirée de ce template  » leaflet timeline » disponible sur github ici. Pour rajouter des graphiques je vous renvoie vers les appi Google chart

Voici un exemple d’intégration

<script type="text/javascript">
  function drawVisualization() {
    //Raw data
    var data = google.visualization.arrayToDataTable([
      ['Années', 'Bureau de seconde main', 'Bureaux neufs',  'Bureaux neufs (Comptes propres et clé en main)'],
      ['2007',  39333,          48511,           0],
     ----
    ]);

    var options = {
      title : 'Transactions neuves et de seconde main de 2007 à 2015 en milliers en m²',
      vAxis: {title: "Surface"},
      //legend: { position: "bottom" },
      //Horizontal axis text vertical
      hAxis: {title: "Années", slantedText:false, slantedTextAngle:10},
      seriesType: "bars",
      //series: {4: {type: "line"}},
      isStacked: true,
       colors:['#3DC2D5','#EC3858','#E83043','#C62939'],
      backgroundColor: { fill:'transparent' },
      fontSize: 10,
    legend: {
    position: 'none',
    alignment: 'start',
    textStyle: {
     color: 'black',
     fontName: 'times',
     fontSize: 10
    } },
   tooltip: {
    isHtml: true,
    showColorCode: true,
    isStacked: true
   }
    };

   // var chart = new google.visualization.BarChart(document.getElementById('chart_div'));
     var chart = new google.visualization.ColumnChart(document.getElementById('chart_div'));

    
    chart.draw(data, options);
  }
  google.setOnLoadCallback(drawVisualization);
</script>

 

#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