OK, así que un segundo intento de responder a su pregunta con una solución GDAL pura.
En primer lugar, GDAL (Biblioteca de abstracción de datos geoespaciales) era originalmente solo una biblioteca para trabajar con datos geoespaciales ráster, mientras que la biblioteca separada de OGR estaba destinada a trabajar con datos vectoriales. Sin embargo, las dos bibliotecas ahora están parcialmente fusionadas, y generalmente se descargan e instalan juntas bajo el nombre combinado de GDAL. Entonces la solución realmente cae bajo OGR. Tienes esto en tu código inicial, así que supongo que lo sabías, pero es una distinción importante para recordar al buscar consejos y sugerencias.
Para leer datos de una capa vectorial, su 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()
Necesitamos crear una nueva característica antes de poder escribirla en un archivo de forma (o cualquier otro conjunto de datos vectoriales). Para crear una nueva entidad, primero necesitamos: - Una geometría - Una definición de entidad, que probablemente incluirá definiciones de campo. Use el constructor Geometry ogr.Geometry () para crear un objeto Geometry vacío. Defina cuál es la geometría de una manera 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()
He tenido suerte leyendo y escribiendo en capas. Específicamente, tengo un código que leerá una capa de archivo de forma que contiene polilíneas y generará la geometría de cada entidad en archivos de texto (utilizados como entrada para un modelo antiguo).
Parece que podría ser útil obtener cada una de las características de sus capas.
Escribir en otra capa no debería ser demasiado complejo desde aquí. Algo como esto debería funcionar en teoría:
Desde aquí, debería poder obtener datos de cada entidad y escribir nuevas entidades en una nueva capa.
Dan
fuente