bonsoir,
j'ai un toc, ce sont les batiments chevauchants. Voila c'est dit!
ci-joint un blabla sur le sujet (il vaut que ce qu'il vaut)
et un script python a tester sur les fichiers cadastre.openstreetmap.fr
Les dépendances sont
+ Osmsax d'osmose :
https://gitorious.org/osmose/backend/blobs/master/analyser_sax.py
shapely (python-shapely)
rtree : http://pypi.python.org/pypi/Rtree
rtree necessite libspatialindex http://libspatialindex.github.com/
Merci de vos remarques
didier
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# building4_qadastre.py
#
# Copyright 2012 didier <[email protected]>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA 02110-1301, USA.
#
# Programme censé corriger des erreurs de chevauchements dans
# les fichiers généré par qadastre2osm
# (qadastre2osm generate OpenStreetMap files from the cadastre through its PDF export.)
# http://gitorious.org/qadastre/qadastre2osm
#
#
#################################################################################################
import OsmSax, sys, copy, time
from rtree import index
# speedus est integer depuis la version 1.2.10. de shapely
#from shapely import speedups
from shapely.geometry import Point, Polygon
#surface maximum de l'intersection pour une correction
SEUIL = 1.0e-7
def coords(n):
return (n.lat, n.lon)
class Node(object):
def __init__(self, id=None, lon=None, lat=None, tags=None):
self.id = id
if lon != None: self.lon, self.lat = float(lon), float(lat)
if tags:
self.tags = tags
else:
self.tags = {}
self.inWay = set()
self.inRel = set()
class Way(object):
def __init__(self, id=None, nodes=None, tags=None):
self.id = id
if nodes:
self.nodes = nodes
else:
self.nodes = []
if tags:
self.tags = tags
else:
self.tags = {}
class Relation(object):
def __init__(self, id, members=None, tags=None):
self.id = id
if members:
self.members = members
else:
self.members = []
if tags:
self.tags = tags
else:
self.tags = {}
def __repr__(self):
return "Relation(id=%r, members=%r, tags=%r)" % (self.id, self.members, self.tags)
class CacheDataOsm:
def __init__(self):
self.nods = {}
self.ways = {}
self.rels = {}
def NodeCreate(self, data):
self.nods[data["id"]] = Node(id=data["id"],lon=data["lon"],lat=data["lat"],tags=data["tag"])
def WayCreate(self, data):
self.ways[data["id"]] = Way(id=data["id"],nodes=data["nd"],tags=data["tag"])
def RelationCreate(self, data):
self.rels[data["id"]] = Relation(id=data["id"],tags=data["tag"],members=data["member"])
def triliste(xliste):
# xliste est une liste (liste des nodes d'un way)
# cette fonction renvoie une liste réordonnée
# du plus petit au plus grand
# qui permet
# d'ordonner de façon unique toutes les combinaisons de description d'un way fermé
# la plus petite valeur de la liste
miref=min(xliste)
# la plus grande valeur de la liste
maref=max(xliste)
# premier tri de la liste
# position de l'element miref dans la liste
pmiref=xliste.index(miref)
# réordonnage de la premiere liste
newl1=xliste[pmiref:len(xliste)]+xliste[1:pmiref+1]
# ecart de position des élément minref et maxref dans la premiere liste
opti1=newl1.index(maref)-newl1.index(miref)
# creation liste temporaire (pour inversion de l'ordre)
tmpl=xliste[pmiref:len(xliste)]+xliste[1:pmiref+1]
# inversion de l'ordre des éléments de la liste
tmpl.reverse()
# position de l'element miref dans la liste temporaire
tpmiref=tmpl.index(miref)
# deuxieme tri de la liste (tri de la liste temporaire)
# deuxieme liste réordonnée
newl2=tmpl[tpmiref:len(tmpl)]+tmpl[1:tpmiref+1]
# critere de description (a^2 + b^3 + ...+ a^n)
opti2=0
cpt=1
for k2 in newl2 :
cpt += 1
opti2 += k2^cpt
# il est peut probable que opti1=opti2
# donc quelque soit la facon dont le way est decrit
# cette fonction retournera toujours la meme liste
if opti1>opti2 :
return newl2
else :
return newl1
def listetostring(xliste) :
stringx=''
for k in xliste :
stringx=stringx+str(k)+";"
return stringx
###########################################################################
avant = time.clock()
# fichier a analyser / corriger en premier argument de la ligne de commande
data = OsmSax.OsmSaxReader(sys.argv[1])
#fout = sys.argv[2]
#data = OsmSax.OsmSaxReader("Y0691-YERRES-houses.osm")
#data = OsmSax.OsmSaxReader("D0229-EVREUX-houses.osm")
fout = 'chevauche-corrige.osm'
dataosm = CacheDataOsm()
print 'Parse du fichier...'
data.CopyTo(dataosm)
print "nombre de node", len(dataosm.nods)
print "nombre de way ", len(dataosm.ways)
# Initialisation des variables
lstCoordId = {}
lstIdCoord = {}
lstPtFinal = {}
lstBounds = {}
listeway={}
c1=[]
c2=[]
nxc1=[]
nxc2=[]
print 'Indexation node...'
# pour retrouver l'id d'un tuple de coordonnée
# pour retrouver le tuple de coordonnée d'un id
# ceci pour les nodes existants ou créés
minid=-1
for xy in dataosm.nods.keys():
# pseudo index
lstCoordId[coords(dataosm.nods[xy])]=dataosm.nods[xy].id
lstIdCoord[dataosm.nods[xy].id]=coords(dataosm.nods[xy])
# pour connaitre l'id minimum des node de ce fichier
if (dataosm.nods[xy].id < minid) :
minid=dataosm.nods[xy].id
# raw_input("Appuyer sur une touche ")
#initialisation du compteur de nouveau node
if minid <0:
nxnx=-minid+1
else :
nxnx=1
print "premier numero des nodes créés:", -nxnx
print 'Indexation way....'
idx = index.Index()
# balayage de tous les ways
cptidx=0
for wx in dataosm.ways.values() :
cptidx+=1
minlat=99999
minlon=99999
maxlat=-99999
maxlon=-99999
for pt in wx.nodes :
if dataosm.nods[pt].lat< minlat : minlat=dataosm.nods[pt].lat
if dataosm.nods[pt].lon< minlon : minlon=dataosm.nods[pt].lon
if dataosm.nods[pt].lat> maxlat : maxlat=dataosm.nods[pt].lat
if dataosm.nods[pt].lon> maxlon : maxlon=dataosm.nods[pt].lon
lstBounds[wx.id]=(minlon,minlat,maxlon,maxlat)
idx.insert(cptidx,lstBounds[wx.id],wx.id)
print 'Correction des chevauchements... patience !'
methode="PointDedans"
for w1 in dataosm.ways.values() :
# coordonnées et id des noeuds
del(c1[:])
#récupération des coordonnées des points du way
for a in w1.nodes:
#c1.append(coords(dataosm.nods[xy]))
c1.append(lstIdCoord[a])
# construction surface
s1= Polygon(c1)
# vérification sur les buildings, pas les inner sans attributs => a corriger manuellement
if not ("building" in w1.tags) : continue
# verification si le way est fermé (le premier et dernier ref du way sont le même node)
if not (w1.nodes[0] == w1.nodes[len(w1.nodes)-1]) : continue
#verification que la forme est géometriquement valide
if (len(w1.nodes)<4) : continue
if s1.is_valid and s1.geom_type=="Polygon":
# liste de tous les way contenus dans le bbox de w1
listew2 = list(idx.intersection(lstBounds[w1.id], objects="raw"))
for widx in listew2:
w2=dataosm.ways[widx]
# pour ne pas traiter son propre chevauchement
if (w1.id <> w2.id) :
if len(w2.nodes)<4 : continue
# verification si le way est fermé (le premier et dernier ref du way sont le même node)
if not (w2.nodes[0] == w2.nodes[len(w2.nodes)-1]) : continue
if not ("building" in w2.tags) : continue
# coordonnées et id des noeuds
del(c2[:])
#récupération des coordonnées des points du way
for b in w2.nodes :
#c2.append(coords(dataosm.nods[xy]))
c2.append(lstIdCoord[b])
# construction surface
s2= Polygon(c2)
#verification que la forme est géometriquement valide
if not s2.is_valid : continue
if s2.geom_type=="Polygon" :
# les 2 formes ont une intersection
if (s1.intersects(s2)==True) :
# la surface d'intersection
ix=s1.intersection(s2)
# pour ne pas corriger des chevauchements trop grand
if not (ix.geom_type=="Polygon") : continue
if (ix.area > 0 ) and (ix.area < SEUIL) :
#print w1.id, w2.id, i2.area*10000000
# nombre de node de s1 contenu dans s2 < nombre de node de s2 contenu dans s1
# s1.contains(les Point de s2)>s2.contains(les Point de s1)
if methode=="PointDedans":
# ----------- Methode 1 pour savoir qui modifier
npt1=0
for xy in c2 :
if s1.contains(Point(xy)) :
npt1+=1
npt2=0
for xy in c1 :
if s2.contains(Point(xy)) :
npt2+=1
# le plus grand batiment aura des nodes en plus
# le plus petit batiment aura sa surface réduite, et des nodes seront a supprimer dans josm (validator)
if npt1>npt2:
castraiter=1
else:
castraiter=2
else :
# ----------- Methode 2 pour savoir qui modifier
if s1.area>s2.area:
castraiter=1
else :
castraiter=2
if castraiter==1 :
# le cas 1
# calcul des 3 surfaces dues aux 2 batiments
i1=s1.difference(s2)
i2=s1.intersection(s2)
i3=s2.difference(s1)
if i1.is_valid and i2.is_valid and i3.is_valid:
# nouvelles valeurs pour s1 et s2
nx1=i1.union(i2)
nx2=s2.difference(s1)
else : continue
else :
# le cas 2
i1=s2.difference(s1)
i2=s2.intersection(s1)
i3=s1.difference(s2)
if i1.is_valid and i2.is_valid and i3.is_valid:
# nouvelles valeurs pour s1 et s2
nx1=s1.difference(s2)
nx2=i1.union(i2)
else : continue
# si les nouvelles surfaces sont ok : valides et de type Polygon (pas une ligne et encore moins un multipolygon)
if not nx1.is_valid : continue
if nx2.is_valid :
#jefai1=0
#jefai2=0
#try :
# if nx1.geom_type=="Polygon" :
# jefai1=1
#except:
# #print w1.id, w2.id, ix.area
# jefai1=0
# "un batiment a l'interieur d'un autre . cela cré une exception car nx1 est "empty"=> solution verifier si nx1 est empty
# w1.tags['erreur'] = 'exception'
# w2.tags['erreur'] = 'exception'
#
#try:
# if nx2.geom_type=="Polygon" :
# jefai2=1
#except:
# jefai1=0
# #print w1.id, w2.id, ix.area
# w1.tags['erreur'] = 'exception'
# w2.tags['erreur'] = 'exception'
#
#if jefai1==1 and jefai2==1:
if nx1.is_empty : continue
if nx2.is_empty : continue
if not (nx1.geom_type=="Polygon"): continue
if nx2.geom_type=="Polygon" :
# retrouver l'id des noeuds des surfaces générées
#---------------------- --- surface 1 -
# liste des coordonnées de la surface1
del(nxc1[:])
for c in nx1.exterior.coords :
try :
# exception si xy n'existe pas dans lstCoordId => Nouveau node
lstPtFinal[c]=lstCoordId[c]
except :
nxnx+=1
lstPtFinal[c]=-nxnx
# ajout des nouveaux noeuds
# mise a jour des dictionnaires coordonnees:point et point:coordonnees
lstCoordId[c]=-nxnx
lstIdCoord[-nxnx]=c
# liste des nodes du nouveau way
for e in nx1.exterior.coords :
nxc1.append(lstCoordId[e])
#--------------------------- surface 2
del(nxc2[:])
for f in nx2.exterior.coords :
try :
# exception si xy n'existe pas dans lstCoordId => Nouveau node
lstPtFinal[f]=lstCoordId[f]
except :
nxnx+=1
lstPtFinal[f]=-nxnx
# ajout des nouveaux noeuds
# mise a jour des dictionnaires coordonnees:point et point:coordonnees
lstCoordId[f]=-nxnx
lstIdCoord[-nxnx]=f
# liste des nodes du nouveau way
for h in nx2.exterior.coords :
nxc2.append(lstCoordId[h])
# ------------------------------ mise a jour des donnees des ways
# modifier les nd des ways
s1= Polygon(nx1.exterior.coords)
w1.nodes=copy.deepcopy(nxc1)
w2.nodes=copy.deepcopy(nxc2)
# ajout tag pour vérification: a supprimer dans josm ou dans le script...
w1.tags['barrier'] = 'wall'
w2.tags['barrier'] = 'wall'
# Corrections du fichier des ways pour optimiser le validator JOSM
doublon={}
nberr1=0
nberr2=0
for w in dataosm.ways.values():
# verification si le way est fermé (le premier et dernier ref du way sont le même node)
if w.nodes[0] == w.nodes[len(w.nodes)-1] :
if len(w.nodes)<4 :
# ajout d'un tag sans attribu pour générer une erreur par le validator de josm
w.tags['fixme']=''
w.tags['note']='way ferme avec moins de 2 noeuds , ce n''est pas une surface'
nberr1 += 1
else :
# nouvelle liste réordonnée
newl=triliste(w.nodes)
# affectation a la liste existante des nouvelles valeurs
w.nodes=copy.deepcopy(newl)
#liste des nodes separés par ;
tmps=listetostring(w.nodes)
# verification si la liste des nodes du way est déja dans la liste => doublon
if tmps in doublon :
# ajout du tag pour identifier le doublon dans josm
# attention pour les doubles inner dans les relations
# par exemple un buiding=yes contenu dans un building= yes et wall=no
nberr2 += 1
#w.tags['fixme']='doublon'
else :
# sinon on ajoute l'element
doublon[tmps]=1
print "Il y a %s erreur(s) et %s doublon(s)." % (nberr1,nberr2)
print 'Ecriture...'
out = OsmSax.OsmSaxWriter(fout, "UTF-8")
out.startDocument()
out.startElement("osm", {'version':'0.6','generator':'ScriptDidier'})
# ajout des nouveaux node
cpt=0
for clef, valeur in lstCoordId.items() :
if valeur <= minid :
# il reste une faille: les noeuds non utilisé dans les ways...
# il suffit de lancer le validator de josm et de virer les erreurs
out.NodeCreate({'id':valeur,'lat':clef[0],'lon':clef[1],'tag':{'fixme':'ajout par script'}})
cpt+=1
print "nombre de node ajoutés: ", cpt
# réécriture des données existantes - existantes modifiées
for n in dataosm.nods.values():
if n.id> minid :
out.NodeCreate({'id':n.id,'lon':n.lon,'lat':n.lat,'tag':n.tags})
else : cpt+=1
print "nombre de node supprimés: ", cpt
for w in dataosm.ways.values():
out.WayCreate({'id':w.id,'nd':w.nodes,'tag':w.tags})
for r in dataosm.rels.values():
out.RelationCreate({'id':r.id,'member':r.members,'tag':r.tags})
out.endElement("osm")
print 'Temps d''execution : ',time.clock() - avant
Préambule:
je n'incrimine personne ou auqu'un outil.
J'énumère des points constatés récurrents par analyse quantitative,
correction ou visualisation
Les requetes effectuées sont un extract france de geofabrik.
Le schéma utilisé est celui d'osm2pgsl- ce qui fausser les chiffres
(relation / batiment , layer distinct, ...)
Données utilisées mars 2012:
france.osm = 43.6 go, dont 31.6 go pour les batiments ( 11 pour les way
et 20 pour les node).
soit 25 224 467 batiments dont 36 mille avec des relations
Bâtiments chevauchants:
Pourquoi
- Ils existent dans les données issues des scripts qui permettent
l'import semi-automatique des bâtiments
+ la précision du cadastre ?, la densité de batiments
+ le type de ligne délimitant les parcelles
+ toutes les erreurs détéctées par osmose après upload ne le
sont pas par le validator de josm avant l'upload
+ plusieurs relations pour le meme outer mais chacune ayant un
inner différent (probleme évoqué sur la liste talk-fr, a prioris corrigé?)
+ ways en doubles (peu nombreux mais dont la correction est
facile)
- Apres l'upload des données, des bâtiments se chevauchent
+ il y a des superposition avec des batiment qui existaient
auparavent
++ dans la meme commune
++ au limites des communes
+ fausse manipulation/problème lors de l'upload qui font que
les batiments sont en plusieurs exemplaire ?
+ ...
Analyse quantitative
Pour la france, le taux d'erreur est de 0.9% (nombre total de
batiments/ nombre d'erreur détecté par osmose)
18 millions de batiments et plus de 200 mille chevauchements. (7
millions sans chevauchement)
Surface de l'intersection
82% fait moins de 2 m2
4% est compris entre 2 et 10 m2
le reste varie jusqu'a 2884m2.
Les plus grosses valeurs concernent zone commerciale, hopital, batiment
industriel
(les batiments avec de grande surface peuvent être a l'origine de
plusieurs erreurs de chevauchement détectées par osmose)
Nombre de Points d'intersection
1 point : 103009 => overlaps 1 node ?
2 points : 37767
3 points : 34816
4 points : 5434
5 points : 14595
6 points : 2750
7 points : 14595
jusqu'a 184 points d'interection !
Solutions:
+ Mettre en avant l'api openstreetmap.fr de sly: auparavant, il était
fastidieux de retélécharger un upload du quadastre , ce qui a pu décourager
certains de la correction post upload
+ voir avec l'equipe d'osmose s'il est possible de distinguer la
"gravité" du chevauchement (critére de surface ou rapport entre la surface de
l'intersection et la plus petite surface des batiments concernés
+ des scripts de correction qui se situerait entre Qadastre2osm et Josm.
+ ...
_______________________________________________
dev-fr mailing list
[email protected]
http://lists.openstreetmap.org/listinfo/dev-fr
