Python GDAL – Splitting Shapefile Per Feature

gdalosgeopythonrasterizationsplit-by-attribute

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

Best Answer

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