Dividindo shapefile por recurso em Python usando GDAL?

é possível dividir um shapefile por recurso em python? (o melhor seria uma solução em que eu possa salvar temporariamente os objetos vetoriais resultantes na memória e no disco).

O motivo: quero usar a função gdal rasterizeLayer com vários subconjuntos diferentes de shapefile. A função requer um objeto osgeo.ogr.Layer.

mkay, tentei um pouco e poderia funcionar da seguinte maneira. Você pode obter a geometria dos objetos da camada gdal por recurso da seguinte maneira.

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

Agora eu só preciso saber como criar um objeto osgeo.ogr.layer com base nessa geometria.

Para esclarecimentos. Eu preciso de uma função no código ogr / gdal simples! Isso parece ser de algum interesse para outras pessoas também e eu ainda quero uma solução sem módulos secundários (embora qualquer solução proveniente daqui seja usada em um plugin qgis disponível gratuitamente).

OK, então uma segunda tentativa de responder sua pergunta com uma solução GDAL pura.

Em primeiro lugar, o GDAL (Geospatial Data Abstraction Library) era originalmente apenas uma biblioteca para trabalhar com dados geoespaciais rasterizados, enquanto a biblioteca OGR separada destinava-se a trabalhar com dados vetoriais. No entanto, as duas bibliotecas agora estão parcialmente mescladas e geralmente são baixadas e instaladas juntas com o nome combinado de GDAL. Portanto, a solução realmente se enquadra no OGR. Você tem isso no seu código inicial, então acho que sabia disso, mas é uma distinção importante a ser lembrada ao procurar dicas e sugestões.

Para ler dados de uma camada vetorial, seu código inicial é bom:

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

Precisamos criar um novo recurso antes de poder gravá-lo em um shapefile (ou qualquer outro conjunto de dados vetoriais). Para criar um novo recurso, primeiro precisamos: - Uma geometria - Uma definição de recurso, que provavelmente incluirá definições de campo Use o construtor Geometry ogr.Geometry () para criar um objeto Geometry vazio. Defina qual é a geometria de maneira diferente para cada tipo (ponto, linha, polígono etc.). Então, por exemplo:

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

ou

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

Para uma definição de campo

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

Agora você pode criar sua camada vetorial. Nesse caso, um polígono quadrado:

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

Tive sorte lendo e escrevendo em camadas. Especificamente, eu tenho um código que lerá uma camada de shapefile contendo polilinhas e produzirá a geometria de cada recurso em arquivos de texto (usados ​​como entrada para um modelo antigo).

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

Parece que pode ser útil obter cada um dos recursos de suas camadas.

Gravar em outra camada não deve ser muito complexo a partir daqui. Algo assim deve funcionar em teoria:

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

A partir daqui, você poderá obter dados de cada recurso e gravar novos recursos em uma nova camada.

Dan