12

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

1

4 Answers 4

9

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

0
1
+50

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!

0
-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

-1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.