Czy podzielić funkcję podczas przecinania się z funkcją innej warstwy za pomocą PyQGIS / Python?

12

Mam warstwę buforową (zielony wielokąt), którą chcę podzielić na dwa wielokąty, gdy tylko przekroczy barierę (niebieska linia). Próbowałem użyć metody „splitGeometry”, ale po prostu nie mogę jej uruchomić. Mój kod do tej pory jest następujący:

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) 

Który zwraca 1 dla wyniku (błąd) i pustą listę dla nowychGeometries. Każda pomoc jest mile widziana.

wprowadź opis zdjęcia tutaj

Alex
źródło
1
Może ten tutaj ci pomoże: gis.stackexchange.com/questions/66543/erase-method-using-ogr
Michalis Avraam

Odpowiedzi:

7

Możesz użyć do tego reshapeGeometryfunkcji QgsGeometryobiektu, który przecina wielokąt wzdłuż jego przecięcia z linią.

Poniższe elementy przecinają wielokąty buforowe z liniami i dodają funkcje podzielonego wielokąta do warstwy pamięci (składnia 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
źródło
1
Działa to doskonale. Najpierw wypróbowałem inne rozwiązanie i zadziałało, więc dałem za to nagrodę, jeszcze zanim przeczytałem ur. To rozwiązanie jest absolutnie idealne i lepiej pasuje do mojego skryptu. przepraszam za to: /
Alex
Hehe, nie ma problemu! Cieszę się, że to pomaga!
Jake,
Głosuję za odpowiedzią, ponieważ działa idealnie, a moja jest jedynie przybliżeniem. @PeyMan Dzięki za nagrodę, ale nie było odpowiedzi oprócz mojej, kiedy wartość nagrody się skończyła. Lepsze rozwiązania są zawsze mile widziane.
Antonio Falciano,
czy jest jakiś sposób na podzielenie całego wielokąta warstwy speicifc?
Muhammad Faizan Khan
Mam jedną warstwę i jest wiele wielokątów, chcę je podzielić przez kodowanie
Muhammad Faizan Khan
2

Dobre przybliżenie za pomocą GDAL> = 1.10.0 skompilowanego z SQLite i SpatiaLite polega na zawijaniu twoich warstw (np. Poligon.shp i line.shp ) w plik OGR VRT (np. 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>

aby mieć bardzo mały bufor (np. 1 mikron) wokół linii. shp otrzymując warstwę * buffer_line *. Następnie możemy zastosować różnicę symetryczną i różnicę do tych geometrii za pomocą 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

Oczywiście wszystkie te rzeczy są doskonale wykonywalne ze skryptu Python:

os.system("some_command with args")

Mam nadzieję że to pomoże!

Antonio Falciano
źródło
@Jake reshapeGeometry zgłasza wyjątek nieznany błąd. Czy istnieje więc inny sposób sprawdzenia przecięcia wielokąta z polilinią?
user99