Splitting feature when intersecting with feature of another layer using PyQGIS

Question

I have a buffer layer (green polygon) which I want to split to two polygons whenever it crosses a barrier (blue line). I have been trying to use "splitGeometry" method, but I just can't get it to work.

My code so far is this :

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)

Which returns 1 for result( error) and an empty list for newGeometries.

enter image description here

Maybe this one here will help you: http://gis.stackexchange.com/questions/66543/erase-method-using-ogr — Michalis Avraam, Oct 04 '13 at 00:43

Jake · Accepted Answer · 2013-10-09T20:39:31.020

You can use the reshapeGeometry function of the QgsGeometry object for this, which cuts a polygon along its intersection with a line.

The following will intersect the buffer polygons with the lines, and add the split polygon features to a memory layer (QGIS 2.0 syntax):

# 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])

score 1 · Answer 2 · answered Oct 09 '13 at 15:28

A good approximation with GDAL >= 1.10.0 compiled with SQLite and SpatiaLite consists into wrapping your layers (e.g. poligon.shp and line.shp) in an OGR VRT file (e.g. 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>

in order to have a very tiny buffer (e.g. 1 micron) around line.shp obtaining the *buffer_line* layer. Then, we can apply symmetric difference and difference on these geometries using 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

Obviously, all this stuff is perfectly executable from a Python script:

os.system("some_command with args")

Hope this helps!

RamonLopezEscudero · Answer 3 · 2017-09-25T21:10:15.363

-1

The answer of Jake have topological problems with the splitted polygons. I made a solution, hope this can help: https://gist.github.com/RamonLopezEscudero/844c1401f5339143da1b2b5cf7ff27bd

edited Sep 25 '17 at 21:10

answered Sep 22 '17 at 21:54

RamonLopezEscudero

1
1

score -1 · Answer 4 · answered Jul 14 '21 at 07:30

-1

You Can use Intersect tool in the Geoprocessing tab in Arcmap.

answered Jul 14 '21 at 07:30

Ahmad Senousi

139
8

Splitting feature when intersecting with feature of another layer using PyQGIS

4 Answers4

Linked