¿Dividir un shapefile por característica en Python usando GDAL?

Question

¿es posible dividir un shapefile por característica en python? (lo mejor sería una solución en la que pueda guardar temporalmente los objetos vectoriales resultantes en la memoria en lugar de en el disco).

El motivo: quiero utilizar la función gdal rasterizeLayer con varios subconjuntos diferentes de shapefile. La función requiere un objeto osgeo.ogr.Layer.

---

mkay, he probado un poco por ahí y podría funcionar de la siguiente manera. Usted puede obtener la geometría de los objetos de la capa gdal por característica como la siguiente.

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

Ahora sólo necesito saber cómo crear un objeto osgeo.ogr.layer basado en esta geometría.

---

Para aclarar. ¡Necesito una función en código ogr/gdal simple! Esto parece ser de algún interés para otras personas también y todavía quiero una solución sin ningún módulo secundario (aunque cualquier solución que viene de aquí se utilizará en un plugin qgis disponible gratuitamente).

Answer 1

Bien, un segundo intento de responder a su pregunta con una solución GDAL pura.

En primer lugar, GDAL (Geospatial Data Abstraction Library) era originalmente sólo una biblioteca para para trabajar con datos geoespaciales rasterizados, mientras que la biblioteca OGR trabajar con datos vectoriales. Sin embargo, las dos bibliotecas están ahora parcialmente fusionadas, y son y generalmente se descargan e instalan juntas bajo el nombre combinado de GDAL. Así que la solución realmente cae bajo OGR. Usted tiene esto en su código inicial, así que supongo que lo sabía, pero es una distinción importante para recordar cuando se buscan consejos y sugerencias.

Para leer datos de una capa vectorial, tu código inicial está bien:

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

Tenemos que crear una nueva característica antes de poder escribirla en un shapefile (o cualquier otro conjunto de datos vectoriales). Para crear una nueva característica, primero necesitamos - Una geometría - Una definición de rasgo, que probablemente incluirá definiciones de campo Utilizar el constructor Geometry ogr.Geometry() para crear un objeto Geometry vacío. Defina lo que es la geometría de forma diferente para cada tipo (punto, línea, polígono, etc.). Así, por ejemplo:

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

o

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

Para una definición de campo

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

Ahora puedes crear tu capa vectorial. En este caso, un polígono cuadrado:

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

Answer 2

He tenido algo de suerte leyendo de y escribiendo a capas. Específicamente, tengo código que leerá una capa shapefile que contiene polilíneas y la salida de la geometría de cada característica a los archivos de texto (utilizado como entrada para un modelo antiguo).

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

Esto parece que podría ser útil para obtener cada una de las características de sus capas.

Escribir en otra capa no debería ser demasiado complejo a partir de aquí. Algo así debería funcionar en teoría:

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

Desde aquí deberías poder obtener los datos de cada característica y escribir nuevas características en una nueva capa.

Dan