[gdal-dev] Layer.DeleteFeature problems

[Jason Roberts]

Spencer,

Although I have not tried the DeleteFeature method (yet), I do recall that
the ESRI Shapefile driver docs mention this:

"The OGR shapefile driver supports rewriting existing shapes in a shapefile
as well as deleting shapes. Deleted shapes are marked for deletion in the
.dbf file, and then ignored by OGR. To actually remove them permanently
(resulting in renumbering of FIDs) invoke the SQL 'REPACK ' via the
datasource ExecuteSQL() method."

Perhaps that would help.

Best,
Jason

-----Original Message-----
From: gdal-dev-bounces at lists.osgeo.org
[mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Spencer Gardner
Sent: Saturday, March 06, 2010 9:00 AM
To: gdal-dev at lists.osgeo.org
Subject: [gdal-dev] Layer.DeleteFeature problems

I'm writing a script that implements Layer.DeleteFeature on an ESRI
shapefile.  Everything appears to be working - no error messages or
anything, but when the script is finished, no features have been deleted.
On a related note, I have tried to circumvent the problem by simply writing
a new true/false field on the features marked for deletion using
Feature.SetField, but that doesn't seem to be writing anything either.

Running TestCapability on my input layer shows that the feature should be
supported, so the error seems to be a failure to write changes to my file.
Maybe I'm missing something?  I'm not sure where I've gone wrong.  Any
advice is appreciated:

Thanks,
Spencer Gardner

# 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
  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()

  # identify nodes for deletion
  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 at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev