Splitting shapefile per feature in Python using GDAL?

Question

is it possible to split a shapefile per feature in python? (best would be a solution where i can temporarily save the resulting vector objects to the memory instead to disk).

The reason: I want to use the gdal rasterizeLayer function with several different subsets of shapefile. The function requires a osgeo.ogr.Layer object.

---

mkay, i tried a little bit around and it could work as following. You can get the geometry of gdal layer objects per feature as following.

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

Now i just need to know how to create a osgeo.ogr.layer object based on this geometry.

---

For clarification. I need a function in plain ogr/gdal code! This appears to be of some interest to other people as well and i still want a solution without any secondary modules (although any solution coming from here will be used in a free available qgis plugin).

Answer 1

OK so a second attempt to answer your question with a pure GDAL solution.

Firstly, GDAL (Geospatial Data Abstraction Library) was originally just a library for working with raster geo-spatial data, while the separate OGR library was intended to work with vector data. However, the two libraries are now partially merged, and are generally downloaded and installed together under the combined name of GDAL. So the solution really falls under OGR. You have this in your initial code so I guess you knew this but it is an important distinction to remember when searching for tips and hints.

To read data from a vector layer, you're initial code is fine:

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

We need to create a new feature before we can write it to a shapefile (or any other vector data set). To create a new feature, we first need: - A geometry - A feature definition, which will probably include field definitions Use the Geometry constructor ogr.Geometry() to create an empty Geometry object. Define what the geometry is in a different way for each type (point, line, polygon, etc.). So for example:

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

or

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

For a field definition

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

Now you can create your vector layer. In this instance, a square 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()

Answer 2

I've had some luck reading from and writing to layers. Specifically, I have code that will read a shapefile layer containing polylines and output the geometry of each feature to text files (used as input for an old model).

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

This seems like it could be useful to get each of the features from your layers.

Writing to another layer shouldn't be too complex from here. Something like this should work in theory:

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

From here you should be able to get data of each feature from and write new features to a new layer.

Dan