Fractionner une entité lors de l'intersection avec une entité d'une autre couche à l'aide de PyQGIS / Python?

12

J'ai une couche tampon (polygone vert) que je souhaite diviser en deux polygones chaque fois qu'elle franchit une barrière (ligne bleue). J'ai essayé d'utiliser la méthode "splitGeometry", mais je n'arrive pas à la faire fonctionner. Jusqu'à présent, mon code est le suivant:

while ldbuffprovider.nextFeature(feat):
  while barprovider.nextFeature(feat2):
    if feat.geometry().intersects(feat2.geometry()):
        intersection = feat.geometry().intersection(feat2.geometry())
        result, newGeometries, topoTestPoints=feat.geometry().splitGeometry(intersection.asPolyline(),True) 

Ce qui renvoie 1 pour le résultat (erreur) et une liste vide pour newGeometries. Toute aide est grandement appréciée.

entrez la description de l'image ici

Alex
la source
1
Peut-être que celui-ci ici vous aidera: gis.stackexchange.com/questions/66543/erase-method-using-ogr
Michalis Avraam

Réponses:

7

Vous pouvez utiliser la reshapeGeometryfonction de l' QgsGeometryobjet pour cela, qui coupe un polygone le long de son intersection avec une ligne.

Les éléments suivants coupent les polygones tampons avec les lignes et ajoutent les entités de polygone fractionné à une couche mémoire (syntaxe QGIS 2.0):

# Get the dataProvider objects for the layers called 'line' and 'buffer'
linepr = QgsMapLayerRegistry.instance().mapLayersByName('line')[0].dataProvider()
bufferpr = QgsMapLayerRegistry.instance().mapLayersByName('buffer')[0].dataProvider()

# Create a memory layer to store the result
resultl = QgsVectorLayer("Polygon", "result", "memory")
resultpr = resultl.dataProvider()
QgsMapLayerRegistry.instance().addMapLayer(resultl)


for feature in bufferpr.getFeatures():
  # Save the original geometry
  geometry = QgsGeometry.fromPolygon(feature.geometry().asPolygon())
  for line in linepr.getFeatures():
    # Intersect the polygon with the line. If they intersect, the feature will contain one half of the split
    t = feature.geometry().reshapeGeometry(line.geometry().asPolyline())
    if (t==0):
      # Create a new feature to hold the other half of the split
      diff = QgsFeature()
      # Calculate the difference between the original geometry and the first half of the split
      diff.setGeometry( geometry.difference(feature.geometry()))
      # Add the two halves of the split to the memory layer
      resultpr.addFeatures([feature])
      resultpr.addFeatures([diff])

Jake
la source
1
Cela fonctionne avec brio. J'ai d'abord essayé l'autre solution et cela a fonctionné, alors j'ai donné la prime pour cela avant même de lire ur ans. Cette solution est absolument parfaite et convient mieux à mon script. désolé pour cela: /
Alex
Hehe, pas de problème! Heureux que ça aide!
Jake
J'évalue votre réponse car elle fonctionne parfaitement, tandis que la mienne n'est qu'une approximation. @PeyMan Merci pour la prime, mais il n'y avait pas de réponses sauf la mienne quand la valeur de la prime a pris fin. De meilleures solutions sont toujours les bienvenues.
Antonio Falciano
existe-t-il un moyen de diviser tous les polygones d'une couche spécifique?
Muhammad Faizan Khan
J'ai une seule couche et il y a plusieurs polygones, je veux les diviser par codage
Muhammad Faizan Khan
2

Une bonne approximation avec GDAL> = 1.10.0 compilé avec SQLite et SpatiaLite consiste à envelopper vos couches (par exemple poligon.shp et line.shp ) dans un fichier OGR VRT (par exemple couches.vrt ):

<OGRVRTDataSource>
    <OGRVRTlayer name="buffer_line">
        <SrcDataSource>line.shp</SrcDataSource>
        <SrcSQL dialect="sqlite">SELECT ST_Buffer(geometry,0.000001) from line</SrcSQL>
    </OGRVRTlayer>
    <OGRVRTlayer name="polygon">
        <SrcDataSource>polygon.shp</SrcDataSource>
    </OGRVRTlayer>
</OGRVRTDataSource>

afin d'avoir un très petit tampon (par exemple 1 micron) autour de line.shp obtenant la couche * buffer_line *. Ensuite, nous pouvons appliquer la différence symétrique et la différence sur ces géométries en utilisant SpatiaLite:

ogr2ogr splitted_polygons.shp layers.vrt -dialect sqlite -sql "SELECT ST_Difference(ST_SymDifference(g1.geometry,g2.geometry),g2.geometry) FROM polygon AS g1, buffer_line AS g2" -explodecollections

Évidemment, tout cela est parfaitement exécutable à partir d'un script Python:

os.system("some_command with args")

J'espère que cela t'aides!

Antonio Falciano
la source
@Jake reshapeGeometry génère une erreur inconnue d'exception. Existe-t-il un autre moyen de vérifier l'intersection entre le polygone et la polyligne?
user99