Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

EDIT:

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.

SECOND EDIT:

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

share|improve this question
add comment

2 Answers 2

up vote 4 down vote accepted

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()
share|improve this answer
    
thanks Dan! I took a different approach and my QGIS-plugin is already functionable (regarding Spatial queries of raster data). Instead of spliting, i created a subset of the underlying raster. You can find a use case on my blog (tinyurl.com/cy6hs9q). Your answer solves the original question, if i want to split and temporarily save vector features. –  Curlew Dec 19 '12 at 22:08
add comment

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

share|improve this answer
    
hey, thanks. Parts of your code will be useful, if i want to write attributes to my shapes. However as i mentioned above i am only using gdal (specifically the gdal.RasterizeFunction) and unless someone knows how to convert a QgsVectorLayer object to a gdal object and vice versa, this question is still unsolved. –  Curlew Nov 14 '12 at 15:50
    
You didn't mention that you need to do this with QGIS. your initial example appears to be plain vanilla ogr. –  DavidF Nov 14 '12 at 17:44
    
i want to do this in QGIS (i need it as a function for a QGIS plugin), but without relying on the QGIS.core modules. Therefore i need the solution in plain ogr. Dan replied because i mentioned in another post that this code is for a QGIS plugin. –  Curlew Nov 14 '12 at 18:18
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.