Contrôle qualité des limites administratives françaises dans OSM

De cquest le mer 20 nov 2013 - 21:08
Le tracé des limites administratives communales françaises est un travail de fourmis que les contributeurs OSM ont démarré il y a déjà plusieurs années.
La seule source de qualité pour laquelle nous avons une autorisation d'utilisation est le cadastre.
 
En 2011 puis en 2013, l'IGN a bien libéré le GEOFLA™®© puis le Route500™®©, mais ces données sont simplifiées et d'une qualité très inférieure à l'objectif que se sont fixé les fourmis mappeuses.
 
Les 100% étant bientôt atteints, se pose la question de la qualité de ces données et leur vérification. En l'absence d'accès aux données "de référence", il n'est pas possible de faire une étude exhaustive, mais a minima une étude sur des points vérifiables.
 
Le GEOFLA et le Route500 contiennent donc des limites simplifiées, mais cette simplification semble n'avoir était faite que sur les polylignes et pas (ou peu) sur les nœuds d'intersection où les frontières de plusieurs communes se rejoignent.
 
J'ai donc creusé cette idée "pour voir" et voici comment j'ai procédé pour comparer ces nœuds dans OSM et Route500.
 
(si les requêtes SQL ne sont pas votre lecture favorite, vous pouvez zapper à la fin de ce billet)
 
Tout d'abord l'import du Route500 dans une base PostGIS...
 
export SHAPE_ENCODING="ISO-8859-1"
ogr2ogr -t_srs EPSG:900913 -f PostgreSQL PG:dbname=osm LIMITE_ADMINISTRATIVE.SHP -overwrite -nlt GEOMETRY -lco SCHEMA=public
ogr2ogr -t_srs EPSG:900913 -f PostgreSQL PG:dbname=osm NOEUD_COMMUNE.SHP -overwrite -nlt GEOMETRY  -lco SCHEMA=public
ogr2ogr -t_srs EPSG:900913 -f PostgreSQL PG:dbname=osm COMMUNE.SHP -overwrite -nlt GEOMETRY  -lco SCHEMA=public
 
Les géométries sont reprojetées en mercator car c'est une base osm2pgsql destinée à la production de tuiles qui va être utilisée côté OSM.
 
Dans psql on procède à un peu de nettoyage et à la création de quelques index qui seront bien utiles... 
 
ALTER TABLE NOEUD_COMMUNE RENAME TO r500_communes_noeuds;
ALTER TABLE r500_communes_noeuds drop ogc_fid;
ALTER TABLE LIMITE_ADMINISTRATIVE RENAME TO r500_limite_admin;
ALTER TABLE r500_limite_admin drop ogc_fid;
ALTER TABLE LIMITE_ADMINISTRATIVE RENAME TO r500_communes;
ALTER TABLE r500_communes drop ogc_fid;
CREATE TABLE r500_noeuds (pt geometry, nb integer, insee text);
CREATE INDEX on r500_noeuds  using gist(pt);
CREATE INDEX on r500_noeuds (insee);
 
Ensuite il faut extraire ces nœuds d'intersection...
INSERT INTO r500_noeuds select ST_Line_Interpolate_Point(l.wkb_geometry,0) as pt, null, null from r500_limite_admin l;
INSERT INTO r500_noeuds select ST_Line_Interpolate_Point(l.wkb_geometry,1) as pt, null, null from r500_limite_admin l;
 
puis déterminer les communes se rencontrant sur ces nœuds...
INSERT INTO r500_noeuds select pt, count(distinct(insee_comm)), string_agg(distinct(insee_comm),' ' order by insee_comm) as insee from r500_noeuds join r500_communes c on (st_touches(pt,wkb_geometry)) join r500_communes_noeuds n on (n.id_rte500=c.id_rte500) group by pt;
 
et faire un petit nettoyage final:
DELETE FROM r500_noeuds where nb is null;
 
La table ainsi générée contient 72598 nœuds (Route500 2013) avec:
  • le noeud d'intersection
  • le nombre de communes touchant cette intersection: 1 pour les iles, enclaves et exclaves, 2 pour les communes du littoral et frontalières, 3 ou plus pour les autres (il y en a seulement 4 qui sont à l'intersection de5 communes, et 1 à l'intersection de 6 communes)
  • les codes INSEE des communes correspondantes (triés par ordre croissant)
 
Création de la table des nœuds d'intersection OSM
 
La table possède la même structure (noeud, nombre de communes, codes insee).
 
CREATE TABLE osm_noeuds (pt geometry, nb integer, insee text);
CREATE INDEX on osm_noeuds  using gist(pt);
CREATE INDEX on osm_noeuds (insee);
 
Pour accélérer les traitements on va créer une table temporaire ne contenant que les polygones des limites administratives:
 
CREATE TEMP TABLE osm_communes as (select "ref:INSEE" as insee, way FROM planet_osm_polygon WHERE boundary is not null and boundary='administrative' and admin_level='8' and "ref:INSEE" is not null);
CREATE INDEX on osm_communes using gist(way);
 
Pour la suite, je n'ai rien trouvé de très optimisé (pour l'instant) et donc...
 
INSERT INTO osm_noeuds select pt, null, null from (select  st_line_interpolate_point((st_dump(st_linemerge(st_intersection(c.way,l.way)))).geom,0) as pt  from osm_communes l join osm_communes c on (st_touches(c.way,l.way) and c.insee > l.insee)) as lim group by pt;
INSERT INTO osm_noeuds select pt, null, null from (select st_line_interpolate_point((st_dump(st_linemerge(st_intersection(c.way,l.way)))).geom,1) as pt  from osm_communes l join osm_communes c on (st_touches(c.way,l.way) and c.insee > l.insee)) as lim group by pt;
 
Puis on dédoublonne...
 
INSERT INTO osm_noeuds select pt, 0,null from osm_noeuds group by pt;
DELETE FROM osm_noeuds where nb is null;
 
Ne reste plus qu'à reprendre chaque nœud pour trouver les communes qui s'y rejoignent:
 
INSERT INTO osm_noeuds select pt, count(*), string_agg(c.insee,' ' order by c.insee) as insee from osm_noeuds join osm_communes c on (st_touches(pt,way)) group by pt;
DELETE FROM osm_noeuds where nb=0;
 
La structure des tables r500_noeuds et osm_noeuds sont donc identiques et on peut commencer les comparaisons... et sortir par exemple la liste des 100 plus gros écarts :
 
SELECT r.insee, min(st_distance(st_transform(r.pt,2154),st_transform(o.pt,2154))) as d from r500_noeuds r join osm_noeuds o on (r.insee=o.insee) group by r.insee order by d desc limit 100;
 
Résultat
 
Un fichier CSV contenant le résultat de cette dernière requête est disponible et mis à jour chaque nuit sur: http://osm13.openstreetmap.fr/~cquest/openfla/ecarts-osm-route500.csv
 
Il contient donc:
  • l'écart dont la valeur médiane est de 22m (50% des écarts sont inférieurs à 22m)
  • les codes INSEE des communes à l'intersection (séparés par un espace).
Il ne reste plus qu'à analyser plus en détail ces écarts et à vérifier les plus importants... un travail de fourmis, mais ça c'est notre spécialité !!
 
Voici une visualisation colorée de ces écarts produite par une nouvelle couche disponible sur tile.openstreetmap.fr