I'm writing a script with 2 primary layers, a road layer and an intersection
layer. I have designed the script to cycle through the intersections and
remove any that have less than 3 road segments connected to it. I have set a
buffer of 1 foot to account for any misalignment. The script runs just fine
(using almost 24,000 intersection nodes) but doesn't identify any intersections
for removal.
After troubleshooting, I have found that counting the road features after
setting the spatial filter returns an unusually high number of road segments
(usually between 9 and 12) for each intersection so that no intersection
returns less than 3 segments. The same phenomenon occurs if I remove the
1-foot buffer. If I test the layers in ArcView by selecting roads that
intersect a specific node, the result is correct. Is there something in my
code that could be causing this problem? The documentation says the
SetSpatialFilter method "may be inaccurately implemented" - would this be a
symptom of this problem? Is there another way of doing this type of analysis?
Thanks for the help. Please excuse my lack of knowledge as I'm fairly new to
this.
Spencer Gardner
University of Wisconsin-Madison
Department of Urban and Regional Planning
# import modules
import os, sys
try:
from osgeo import ogr
except:
import ogr
startDir = os.getcwd()
tempDir = '/tmp'
# set the working directory
os.chdir(startDir + '/Roads')
# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# open nodes.shp and get the layer
nodesDS = driver.Open('nodesdelete.shp', 0)
if nodesDS is None:
print 'Could not open nodesdelete.shp'
sys.exit(1)
nodesLayer = nodesDS.GetLayer()
# open 2009roads.shp and get the layer
roadsDS = driver.Open('2009roads.shp', 0)
if roadsDS is None:
print 'Could not open 2009roads.shp'
sys.exit(1)
roadsLayer = roadsDS.GetLayer()
print 'found' + str(roadsLayer.GetFeatureCount()) + 'total road features'
# set up the array to store id's to be deleted
deleteNodes = []
# set up the progress bar
count = 1
totalNodes = nodesLayer.GetFeatureCount()
print 'found ' + str(totalNodes) + 'total nodes'
# loop through each node to count the number of connected roads
nodesLayer.ResetReading()
nodesFeature = nodesLayer.GetNextFeature()
while nodesFeature:
# buffer the nodes by 1 foot (to get only the adjacent roads)
nodesGeom = nodesFeature.GetGeometryRef()
bufferGeom = nodesGeom.Buffer(1)
# use bufferGeom as a spatial filter on the node to get all roads
# within 1 foot of the node
roadsLayer.SetSpatialFilter(bufferGeom)
roadCount = roadsLayer.GetFeatureCount()
print str(roadsLayer.GetFeatureCount())
if roadCount < 3:
deleteNodes.append(nodesFeature.GetFID())
# show progress
count += 1
sys.stdout.write('.')
progress = divmod(1000*count, 24000)
if progress[1] == 0:
print str(count) + '/' + str(24000)
# move to the next node
nodesFeature.Destroy()
nodesFeature = nodesLayer.GetNextFeature()
# delete the features
print 'deleting features\n'
for node in deleteNodes:
nodesLayer.DeleteFeature(node)
sys.stdout.write('.')
# close the shapefiles
nodesDS.Destroy()
roadsDS.Destroy()
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
