Membagi shapefile per fitur dalam Python menggunakan GDAL?

apakah mungkin untuk membagi shapefile per fitur dalam python? (terbaik akan menjadi solusi di mana saya dapat menyimpan sementara objek vektor yang dihasilkan ke memori, bukan ke disk).

Alasannya: Saya ingin menggunakan fungsi rasterizeLayer gdal dengan beberapa subset berbeda dari shapefile. Fungsi ini membutuhkan objek osgeo.ogr.Layer.

mkay, saya mencoba sedikit di sekitar dan itu bisa berfungsi sebagai berikut. Anda bisa mendapatkan geometri objek lapisan gdal per fitur sebagai berikut.

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

Sekarang saya hanya perlu tahu cara membuat objek osgeo.ogr.layer berdasarkan geometri ini.

Untuk klarifikasi. Saya membutuhkan fungsi dalam kode ogr / gdal biasa! Ini tampaknya menarik bagi orang lain juga dan saya masih menginginkan solusi tanpa modul sekunder (walaupun solusi apa pun yang berasal dari sini akan digunakan dalam plugin qgis gratis yang tersedia).

OK jadi upaya kedua untuk menjawab pertanyaan Anda dengan solusi GDAL murni.

Pertama, GDAL (Geospatial Data Abstraction Library) pada awalnya hanya perpustakaan untuk bekerja dengan data geo-spasial raster, sedangkan perpustakaan OGR terpisah dimaksudkan untuk bekerja dengan data vektor. Namun, kedua perpustakaan sekarang sebagian digabungkan, dan umumnya diunduh dan diinstal bersama di bawah nama gabungan GDAL. Jadi solusinya benar-benar berada di bawah OGR. Anda memiliki ini dalam kode awal Anda, jadi saya kira Anda tahu ini, tetapi ini merupakan perbedaan penting untuk diingat ketika mencari tips dan petunjuk.

Untuk membaca data dari layer vektor, kode awal Anda baik-baik saja:

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

Kita perlu membuat fitur baru sebelum kita dapat menulisnya ke shapefile (atau kumpulan data vektor lainnya). Untuk membuat fitur baru, pertama-tama kita perlu: - Sebuah geometri - Sebuah definisi fitur, yang mungkin akan menyertakan definisi bidang. Gunakan ogr.Geometry konstruktor Geometri () untuk membuat objek Geometri kosong. Tentukan apa geometri dengan cara yang berbeda untuk setiap jenis (titik, garis, poligon, dll.). Jadi misalnya:

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

atau

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

Untuk definisi bidang

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

Sekarang Anda dapat membuat layer vektor Anda. Dalam hal ini, poligon persegi:

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

Saya sudah beruntung membaca dari dan menulis ke lapisan. Secara khusus, saya memiliki kode yang akan membaca lapisan shapefile yang berisi polyline dan menampilkan geometri setiap fitur ke file teks (digunakan sebagai input untuk model lama).

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

Ini sepertinya bermanfaat untuk mendapatkan setiap fitur dari lapisan Anda.

Menulis ke lapisan lain seharusnya tidak terlalu rumit dari sini. Sesuatu seperti ini seharusnya bekerja secara teori:

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

Dari sini Anda harus bisa mendapatkan data dari setiap fitur dari dan menulis fitur baru ke layer baru.

Dan