Annonce
Pour sa 21ème année, l’association GeoRezo a toujours besoin de vous !
10€ = 1 mois de frais bancaires ; 15€ = 12 mois de nom de domaine ; 30€ = 1 semaine de location des serveurs …
Retrouver nos membres bienfaiteurs
#1 Fri 31 January 2020 15:14
- Olivier Pompier
- Participant occasionnel
- Date d'inscription: 8 Sep 2013
- Messages: 49
QGIS 3/PyQgis: Auto-intersection
Bonjour à la communauté,
Je sollicite votre aide concernant un script qui a pour but d'intersecter les éléments d'une couche et de calculer le nombre de superpositions. J'ai pondu le script suivant qui fait le job:
Code:
while compteur > 0:
outFn = os.path.join(os.path.dirname(path_file),'intersection'+str(compteur)+'.shp')
writer = QgsVectorFileWriter(outFn, 'UTF-8', fields, QgsWkbTypes.Polygon, layer.sourceCrs(), 'ESRI Shapefile')
l1 = list()
for feat in features:
l1.append(feat)
l2 = l1 # la liste 2 est une copie de la liste 1
l3 = [] # liste contenant les surfaces des éléments intersectés
for elt1 in l1:
for elt2 in l2:
if elt1.id() > elt2.id():
if elt1.geometry().intersects(elt2.geometry()):
intersection = elt1.geometry().intersection(elt2.geometry())
elt = QgsFeature(fields)
elt.setAttribute('nb_sup',compteur)
elt.setGeometry(intersection)
if round(elt.geometry().area(),4) not in l3:
writer.addFeature(elt)
l3.append(round(elt.geometry().area(),4))
# iface.addVectorLayer(outFn, '','ogr')
del(writer)
if len(l3) == 0:
compteur = 0
else:
compteur += 1
layer = QgsVectorLayer(outFn, str(compteur), "ogr")
features = layer.getFeatures()Problème, lorsque la couche contient beaucoup d' éléments, le temps de traitement est extrêmement long (des heures et des heures), et généralement, je stoppe l’exécution du script avant.
Auriez vous des suggestions / idées pour un meilleur code ?
Merci, O. Pompier
Hors ligne
#2 Sun 02 February 2020 21:43
- JD
- Moderateur
- Date d'inscription: 8 Aug 2013
- Messages: 726
Re: QGIS 3/PyQgis: Auto-intersection
Bonjour,
je vous invite à utiliser les index spatiaux, cela devrait considérablement améliorer les temps de traitement et faire une seule boucle sur le jeu de donnée.
Un truc dans le genre :
Code:
index = QgsSpatialIndex(features)
for elt1 in l1:
bbox = elt1.geom().boundingBox()
intersects = index.intersects(bbox)
for i in intersects:
if i.id() <> elt1.id():
if elt1.geometry().intersects(i.geometry()):Hors ligne
#3 Thu 06 February 2020 17:28
- Olivier Pompier
- Participant occasionnel
- Date d'inscription: 8 Sep 2013
- Messages: 49
Re: QGIS 3/PyQgis: Auto-intersection
Merci JD pour votre réponse. Effectivement, cela a diminué le temps de traitements significativement.
Le code fonctionnel ci-après:
Code:
index = QgsSpatialIndex(features)
feats = [feat for feat in layer.getFeatures()]
for feat in feats:
for id in index.intersects(feat.geometry().boundingBox()):
if feat.id() > id:
intersection = feat.geometry().intersection(feats[id].geometry())
elt = QgsFeature(fields)
elt.setAttribute('nb_sup',compteur)
elt.setGeometry(intersection)Hors ligne
