That makes sense. I was using FTools so perhaps the GEOS support is not built in? Anyway, I ran it on my Linux box and the analysis part works fine. However, the features that it identifies for deletion don't actually get deleted so the file remains unchanged.
(code included at bottom) The documentation mentions that Layer.DeleteFeature() doesn't work for some drivers, but I don't see any error messages. Does the ESRI Shapefile driver not support deletion? If so, how does one remove features from a shapefile layer? Thanks, Spencer Gardner University of Wisconsin-Madison Department of Urban and Regional Planning ----- Original Message ----- From: Frank Warmerdam <warmer...@pobox.com> Date: Thursday, March 4, 2010 6:11 pm Subject: Re: [gdal-dev] SetSpatialFilter not working? To: Spencer Gardner <sjgard...@wisc.edu> Cc: gdal-dev@lists.osgeo.org > Spencer Gardner wrote: > >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, > > If GEOS is linked into GDAL/OGR the spatial filtering should be accurate. > Otherwise a bounding box comparison is done and you can get significant > numbers of false positives depending on the structure of the data. > > So, my suggestion is to rebuild GDAL with GEOS support. > > Best regards, > -- > ---------------------------------------+-------------------------------------- > I set the clouds in motion - turn up | Frank Warmerdam, warmer...@pobox.com > light and sound - activate the windows | http://pobox.com/~warmerdam > and watch the world go round - Rush | Geospatial Programmer for Rent > # 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', 1) 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 and count < 100: # 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() 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 'identified ' + str(len(deleteNodes)) + ' features to be deleted' 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