¿Dividir una entidad cuando se cruza con una entidad de otra capa usando PyQGIS / Python?

12

Tengo una capa intermedia (polígono verde) que quiero dividir en dos polígonos cada vez que cruza una barrera (línea azul). He estado tratando de usar el método "splitGeometry", pero no puedo hacerlo funcionar. Mi código hasta ahora es este:

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) 

Que devuelve 1 para el resultado (error) y una lista vacía para newGeometries. Cualquier ayuda es muy apreciada.

ingrese la descripción de la imagen aquí

Alex
fuente
1
Quizás este aquí te ayude: gis.stackexchange.com/questions/66543/erase-method-using-ogr
Michalis Avraam

Respuestas:

7

Puede usar la reshapeGeometryfunción del QgsGeometryobjeto para esto, que corta un polígono a lo largo de su intersección con una línea.

Lo siguiente intersectará los polígonos del búfer con las líneas y agregará las características del polígono dividido a una capa de memoria (sintaxis 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
fuente
1
Esto funciona de manera brillante. Primero probé la otra solución y funcionó, así que me dio la recompensa por eso incluso antes de leer sus preguntas. Esta solución es absolutamente perfecta y se adapta mejor a mi script. lo siento por eso: /
Alex
Jeje, no hay problema! Me alegra que ayude!
Jake
Voto su respuesta porque funciona perfectamente, mientras que la mía es solo una aproximación. @PeyMan Gracias por la recompensa, pero no hubo respuestas excepto la mía cuando terminó la recompensa. Las mejores soluciones son siempre bienvenidas.
Antonio Falciano
¿Hay alguna manera de dividir todo el polígono de una capa específica?
Muhammad Faizan Khan
Tengo una sola capa y hay varios polígonos, quiero dividirlos a través de la codificación
Muhammad Faizan Khan
2

Una buena aproximación con GDAL> = 1.10.0 compilado con SQLite y SpatiaLite consiste en envolver sus capas (por ejemplo, poligon.shp y line.shp ) en un archivo OGR VRT (por ejemplo, layers.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>

para tener un búfer muy pequeño (por ejemplo, 1 micrón) alrededor de line.shp obteniendo la capa * buffer_line *. Luego, podemos aplicar la diferencia simétrica y la diferencia en estas geometrías usando 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

Obviamente, todo esto es perfectamente ejecutable desde un script Python:

os.system("some_command with args")

¡Espero que esto ayude!

Antonio Falciano
fuente
@Jake reshapeGeometry está lanzando un error desconocido de excepción. Entonces, ¿hay alguna otra forma de verificar la intersección entre el polígono y la polilínea?
user99