I have two shapefiles that intersect: one containing a line, and the other containing polygons. I'm trying to use OGR and Python to create a buffer around the line and select polygons within it, then write those features to a shapefile.
##import modules
import ogr
import gdalconst
import sys
import os
##Define input files
power = (r'C:\Users\Megan\Desktop\powerlines\PowerLine.shp')
parcel = (r'C:\Users\Megan\Desktop\parcels\Parcels.shp')
outName = ("BufferedParcels.shp")
##Define buffer width
bufferWidth = 250
##Assign driver
driver = ogr.GetDriverByName('ESRI Shapefile')
##Open files
powerData = driver.Open(power, gdalconst.GA_ReadOnly)
parcelData = driver.Open(parcel, gdalconst.GA_ReadOnly)
##Verify files were opened correctly, exit if not
if powerData is None:
print('Could not open ' + power)
sys.exit(1)
if parcelData is None:
print('Could not open ' + parcel)
sys.exit(1)
inputDS = driver.Open(power)
if inputDS is None:
print('Could not open', power)
sys.exit(1)
##Get feature layer and geometry from powerlines
powerLayer = powerData.GetLayer(0)
powerFeature = powerLayer.GetNextFeature()
powerGeom = powerFeature.GetGeometryRef()
##Get feature layer and count from parcels
parcelLayer = parcelData.GetLayer(0)
parcelFeatureCo = parcelLayer.GetFeatureCount()
##Re-write output file if it already exists
inputLayer = inputDS.GetLayer()
if os.path.exists(outName):
os.remove(outName)
##Apply buffer to powerline geometry, give affordance that it worked
powerBuffer = powerGeom.Buffer(bufferWidth)
print('Created 250 ft buffer around',power)
##Attempt to get output driver, exit in case of failure
try:
outputDS = driver.CreateDataSource(outName)
except:
print('Could not create',outName)
sys.exit(1)
##Define new layer, exit in case of failure
newLayer = outputDS.CreateLayer('BufferedParcels',geom_type = ogr.wkbPolygon)
if newLayer is None:
print('Could not create layer for buffer in output DS')
sys.exit(1)
##Get new layer definition
newLayerDef = newLayer.GetLayerDefn()
##Give affordance that program is running
##Set feature ID
##Loop through features in parcelLayer
print ('Writing the following parcels to shapefile:')
featureID = 0
for i in parcelLayer:
##parcelFeatures = parcelLayer.GetFeature(i)
parcelGeom = i.GetGeometryRef()
bufferContains = powerBuffer.Contains(parcelGeom)
names = i.GetField('SITUSADDR')
area = i.GetField('AREA')
relevantVals = (names, area)
fmt = '%s %s'
##Determine if features are within buffer
##Print stats if 'True'
if bufferContains == True:
print (fmt % relevantVals)
try:
newFeature = ogr.Feature(newLayerDef)
newFeature.SetGeometry(poly)
newFeature.SetFID(featureID)
newLayer.CreateFeature(newFeature)
except:
##print('Error printing shapefile.')
newFeature.Destroy()
##parcelLayer.Destroy()
featureID += 1
i.Destroy()
##Give affordance that program worked
print('Wrote', featureID, 'parcels to shapefile.')
inputDS.Destroy()
outputDS.Destroy()
The code above successfully writes a list of the polygons I'm looking for and creates a shapefile. However, the file has an empty attribute table and doesn't show any geometry.
Can anyone see where I'm going wrong?
1 Answer 1
It seems that you missed the correct variable in the statement in the code above:
newFeature.SetGeometry(poly)
which should be:
newFeature.SetGeometry(parcelGeom)
And as @klewis suggested, adding attributes using SetField should help populate the attribute table.
Additionally, you can also try putting the datasource out of scope in a different manner using outputDS = None instead of using Destroy().
NameError: name 'outputDS' is not defined...?