Teilen eines Shapefiles pro Feature in Python mit GDAL?

14

Ist es möglich, ein Shapefile pro Feature in Python zu teilen? (Am besten wäre eine Lösung, bei der ich die resultierenden Vektorobjekte vorübergehend im Speicher anstatt auf der Festplatte speichern kann).

Der Grund: Ich möchte die Funktion gdal rasterizeLayer mit mehreren unterschiedlichen Teilmengen von Shapefile verwenden. Die Funktion benötigt ein osgeo.ogr.Layer-Objekt.


mkay, ich habe ein bisschen rumprobiert und es könnte so funktionieren. Sie können die Geometrie von GDAL-Layer-Objekten pro Feature wie folgt abrufen.

#  Load shape into gdal
shapefile=str(vectorPath)
layer_source = ogr.Open(shapefile)
lyr = layer_source.GetLayer(0)
for i in range(0,lyr.GetFeatureCount()):
     feat = lyr.GetFeature(i)
     ge = feat.geometry()

Jetzt muss ich nur noch wissen, wie man ein osgeo.ogr.layer-Objekt basierend auf dieser Geometrie erstellt.


Zur Klarstellung. Ich brauche eine Funktion in normalem OGR / GDAL-Code! Dies scheint auch für andere von Interesse zu sein und ich möchte immer noch eine Lösung ohne sekundäre Module (obwohl jede von hier kommende Lösung in einem frei verfügbaren QGIS-Plugin verwendet wird).

Brachvogel
quelle

Antworten:

7

OK, also ein zweiter Versuch, Ihre Frage mit einer reinen GDAL-Lösung zu beantworten.

Erstens war GDAL (Geospatial Data Abstraction Library) ursprünglich nur eine Bibliothek für die Arbeit mit georäumlichen Raster-Daten, während die separate OGR-Bibliothek für die Arbeit mit Vektordaten gedacht war. Die beiden Bibliotheken werden jetzt jedoch teilweise zusammengeführt und in der Regel zusammen unter dem kombinierten Namen GDAL heruntergeladen und installiert. Die Lösung fällt also wirklich unter OGR. Sie haben dies in Ihrem ursprünglichen Code, also denke ich, dass Sie das gewusst haben, aber es ist ein wichtiger Unterschied, den Sie beachten müssen, wenn Sie nach Tipps und Hinweisen suchen.

Zum Lesen von Daten aus einer Vektorebene reicht der ursprüngliche Code aus:

from osgeo import ogr
shapefile = ogr.Open(shapefile)
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    print i, name, geometry.GetGeometryName()

Wir müssen ein neues Feature erstellen, bevor wir es in ein Shapefile (oder einen anderen Vektordatensatz) schreiben können. Um ein neues Feature zu erstellen, benötigen wir zunächst: - Eine Geometrie - Eine Feature-Definition, die wahrscheinlich Felddefinitionen enthält. Verwenden Sie den Geometriekonstruktor ogr.Geometry (), um ein leeres Geometrieobjekt zu erstellen. Definieren Sie die Geometrie für jeden Typ unterschiedlich (Punkt, Linie, Polygon usw.). Also zum Beispiel:

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,20)

oder

line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10,10)
line.AddPoint(20,20)
line.SetPoint(0,30,30) #(10,10) -> (30,30)

Für eine Felddefinition

fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)

Jetzt können Sie Ihre Vektorebene erstellen. In diesem Fall ein quadratisches Polygon:

#create simple square polygon shapefile:
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')

datasource = driver.CreateDataSource('YOUR_PATH')
layer = datasource.CreateLayer('layerName',geom_type=ogr.wkbPolygon)

#create polygon object:
myRing = ogr.Geometry(type=ogr.wkbLinearRing)
myRing.AddPoint(0.0, 0.0)     #LowerLeft
myRing.AddPoint(0.0, 10.0)    #UpperLeft
myRing.AddPoint(10.0, 10.0)   #UpperRight
myRing.AddPoint(10.0, 0.0)    #Lower Right
myRing.AddPoint(0.0, 0.0)     #close ring
myPoly = ogr.Geometry(type=ogr.wkbPolygon)
myPoly.AddGeometry(myRing)
print ('Polygon area =',myPoly.GetArea())  #returns correct area of 100.0

#create feature object with point geometry type from layer object:
feature = ogr.Feature( layer.GetLayerDefn())
feature.SetGeometry(myPoly)
layer.CreateFeature(feature)

#flush memory - very important
feature.Destroy()
datasource.Destroy()
Dan
quelle
Danke Dan! Ich habe einen anderen Ansatz gewählt und mein QGIS-Plugin ist bereits funktionsfähig (in Bezug auf räumliche Abfragen von Rasterdaten). Anstatt zu splitten, habe ich eine Teilmenge des zugrunde liegenden Rasters erstellt. Einen Anwendungsfall finden Sie in meinem Blog ( tinyurl.com/cy6hs9q ). Ihre Antwort löst die ursprüngliche Frage, ob ich Vektormerkmale teilen und vorübergehend speichern möchte.
Curlew
5

Ich hatte etwas Glück beim Lesen und Schreiben in Schichten. Insbesondere habe ich Code, der eine Shapefile-Ebene mit Polylinien liest und die Geometrie jedes Features in Textdateien ausgibt (die als Eingabe für ein altes Modell verwendet werden).

name     = layer.name()
provider = layer.dataProvider()
feat     = QgsFeature()

# Now we can loop through all the defined features
while provider.nextFeature(feat):

    # Get layer attributes               
    attrs = feat.attributeMap()
    for (k,attr) in attrs.iteritems():
        if k == 0:
            attrOne = attr.toString()
        elif k == 1:
            attrTwo = attr.toString()
        ...

    # Gets the geometry of the feature
    geom = feat.geometry()

    # Get the coordinates of the whole line [or use asPoint()]                    
    line = geom.asPolyline()
        # all points in the line
        for point in line:
            lat = point[0]
            lon = point[1]
            # Add these to a QgsGeometry
            your_Own_QgsGeometry.add...

Dies scheint nützlich zu sein, wenn Sie die einzelnen Features aus Ihren Layern abrufen.

Das Schreiben auf eine andere Ebene sollte von hier aus nicht zu komplex sein. So etwas sollte theoretisch funktionieren:

# New layer name
filename = "myNewLayer.shp"

# define fields for feature attributes
fields   = { 0 : QgsField("attrOne", QVariant.String),
             1 : QgsField("attrTwo", QVariant.String),
             2 : QgsField("...", QVariant.Int) }

# Create coordinate reference system as WGS84
crs    = QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.PostgisCrsId)

# Create the vector layer
writer = QgsVectorFileWriter(filename, "CP1250", fields, QGis.WKBLineString, crs)

# Create some features
feat = QgsFeature()
feat.addAttribute(0, QVariant(runway))
feat.addAttribute(1, QVariant(arriveDepart))
feat.addAttribute(2, QVariant(subTrack))

# Add your geometry
feat.setGeometry(your_Own_QgsGeometry)

# Add the features
writer.addFeature(feat)

# Add it to QGIS project
self.iface.addVectorLayer(filename, "MyLayerName", "ogr")

Von hier aus sollten Sie in der Lage sein, Daten für jedes Feature abzurufen und neue Features in eine neue Ebene zu schreiben.

Dan

Dan
quelle
Hey danke. Teile Ihres Codes werden nützlich sein, wenn ich Attribute zu meinen Formen schreiben möchte. Wie oben erwähnt, verwende ich jedoch nur gdal (speziell die Funktion gdal.RasterizeFunction) und es sei denn, jemand weiß, wie man ein QgsVectorLayer-Objekt in ein gdal-Objekt konvertiert und umgekehrt, ist diese Frage noch ungelöst.
Brachvogel
Sie haben nicht erwähnt, dass Sie dies mit QGIS tun müssen. Ihr erstes Beispiel scheint Vanille Ogr zu sein.
DavidF
Ich möchte dies in QGIS tun (ich brauche es als Funktion für ein QGIS-Plugin), aber ohne auf die QGIS.core-Module angewiesen zu sein. Deshalb brauche ich die Lösung in plain ogr. Dan antwortete, weil ich in einem anderen Beitrag erwähnt habe, dass dieser Code für ein QGIS-Plugin ist.
Brachvogel