Dividir um recurso ao cruzar com um recurso de outra camada usando PyQGIS / Python?

12

Eu tenho uma camada de buffer (polígono verde) que quero dividir em dois polígonos sempre que cruzar uma barreira (linha azul). Eu tenho tentado usar o método "splitGeometry", mas simplesmente não consigo fazê-lo funcionar. Meu código até agora é 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 retorna 1 para resultado (erro) e uma lista vazia para newGeometries. Qualquer ajuda é muito apreciada.

insira a descrição da imagem aqui

Alex
fonte
1
Talvez este aqui irá ajudá-lo: gis.stackexchange.com/questions/66543/erase-method-using-ogr
Michalis Avraam

Respostas:

7

Você pode usar a reshapeGeometryfunção do QgsGeometryobjeto para isso, que corta um polígono ao longo de sua interseção com uma linha.

A seguir, os polígonos do buffer serão cruzados com as linhas e os recursos de polígono dividido serão adicionados a uma camada de memória (sintaxe 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
fonte
1
Isso funciona de maneira brilhante. Eu tentei a outra solução primeiro e ela funcionou, por isso deu a recompensa por isso antes mesmo de ler nossos anseios. Esta solução é absolutamente perfeita e se adapta melhor ao meu script. desculpe por isso: /
Alex
Hehe, não tem problema! Ainda bem que ajuda!
Jake
Votei sua resposta porque ela funciona perfeitamente, enquanto a minha é apenas uma aproximação. @PeyMan Obrigado pela recompensa, mas não havia respostas, exceto as minhas, quando a recompensa acabou. Melhores soluções são sempre bem-vindas.
Antonio Falciano 10/10
existe alguma maneira de dividir todo o polígono de uma camada específica?
Muhammad Faizan Khan
Eu tenho uma única camada e há múltiplos polígono eu quero dividi-los através de codificação
Muhammad Faizan Khan
2

Uma boa aproximação com GDAL> = 1.10.0 compilada com SQLite e SpatiaLite consiste em agrupar suas camadas (por exemplo, poligon.shp e line.shp ) em um arquivo OGR VRT (por exemplo, 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 ter um buffer muito pequeno (por exemplo, 1 mícron) em torno do line.shp, obtendo a camada * buffer_line *. Em seguida, podemos aplicar diferença simétrica e diferença nessas geometrias usando o 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, tudo isso é perfeitamente executável a partir de um script Python:

os.system("some_command with args")

Espero que isto ajude!

Antonio Falciano
fonte
O @Jake reshapeGeometry está lançando um erro desconhecido de exceção.
user99