2

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?

  • Sure thing, added in the code above ^ – Megan Apr 15 at 20:14
  • Still missing something somewhere: NameError: name 'outputDS' is not defined...? – Phil Apr 15 at 20:20
  • Hmm, okay one sec. I define it in the code where the "..." is right now. I'll add the missing blocks (I haven't gotten NameError-ed when I'm running it). – Megan Apr 15 at 20:27
  • Do you need to add, newFeature.SetField('SITUSADDR', names), newFeature.SetField('AREA', area) after SetGeometry? – klewis Apr 15 at 20:33
  • Oh, good thinking @klewis. I just added those lines. It didn't work, I'm still getting an empty attribute table and no geometry. I'm sure that was part of the problem, but I'm stumped on what else it could be! – Megan Apr 15 at 20:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.