Разделить элемент при пересечении с элементом другого слоя, используя PyQGIS / Python?

12

У меня есть буферный слой (зеленый многоугольник), который я хочу разделить на два многоугольника всякий раз, когда он пересекает барьер (синяя линия). Я пытался использовать метод "splitGeometry", но я просто не могу заставить его работать. Мой код пока таков:

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) 

Который возвращает 1 для результата (ошибка) и пустой список для newGeometries. Любая помощь очень ценится.

введите описание изображения здесь

Alex
источник
1
Может быть, этот поможет вам: gis.stackexchange.com/questions/66543/erase-method-using-ogr
Михалис Авраам

Ответы:

7

Для этого вы можете использовать reshapeGeometryфункцию QgsGeometryобъекта, которая разрезает многоугольник вдоль его пересечения с линией.

Следующее будет пересекать буферные полигоны со строками и добавлять элементы разделенного полигона в слой памяти (синтаксис 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])

Джейк
источник
1
Это работает блестяще. Сначала я попробовал другое решение, и оно сработало, поэтому я получил вознаграждение за это еще до того, как прочитал ваш ответ. Это решение идеально подходит и лучше подходит для моего сценария. извините за это: /
Алекс
Хе-хе, нет проблем! Рад, что помогает!
Джейк
Я приветствую ваш ответ, потому что он отлично работает, а мой - только приблизительный. @PeyMan Спасибо за награду, но не было ответов, кроме моего, когда ценность награды закончилась. Лучшие решения всегда приветствуются.
Антонио Фальчано
Есть ли способ разделить весь полигон специального слоя?
Мухаммед Файзан Хан
У меня есть один слой, и есть несколько многоугольников, я хочу разделить их через кодирование
Мухаммед Файзан Хан
2

Хорошее приближение с GDAL> = 1.10.0, скомпилированным с SQLite и SpatiaLite, состоит в том, чтобы обернуть ваши слои (например, poligon.shp и line.shp ) в файл OGR VRT (например, layer.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>

чтобы иметь очень маленький буфер (например, 1 микрон) вокруг line.shp, получая слой * buffer_line *. Затем мы можем применить симметричную разность и разность к этим геометриям, используя 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

Очевидно, что все это прекрасно выполняется из скрипта Python:

os.system("some_command with args")

Надеюсь это поможет!

Антонио Фальчано
источник
@Jake reshapeGeometry вызывает исключительную ошибку исключения. Так есть ли другой способ проверить пересечение между полигоном и полилинией?
user99