How to select polygons by criteria and save them to shp file?

Question

I have file of polygons and want to select few by criteria - area more or less then specific number and save result to shp file.

code gives 0 by passing criteria:

from osgeo import ogr, gdal

shpFile = 'ReBu.shp'

driver = ogr.GetDriverByName('ESRI Shapefile')
ptsDS = driver.Open(shpFile)
ptsLayer = ptsDS.GetLayer()


for i in range(0, ptsLayer.GetFeatureCount()):
    pt = ptsLayer.GetNextFeature()
    print pt.GetField('Area_ha')

myQuery = 'Area_ha = 39.8953'
ptsLayer.SetAttributeFilter(myQuery)

score 2 · Answer 1 · answered Apr 12 '16 at 18:21

It's always risky to attempt equivalence operations on floating-point values. Even if the value displays "39.8953", it's likely to be stored as "39.8952999999999999997" or "39.895300000000000001". The only safe way to compare against floating-point values is perform range-checking:

myQuery = 'Area_ha >= 39.89525 and Area_ha < 39.89535'

Depending on the driver, you might be able to use an absolute value function to change the expression:

myQuery = 'abs(Area_ha - 39.8953) < 0.00005'

though this sort of query would not be able to use an attribute index, and would therefore run much slower.

score 2 · Answer 2 · edited Apr 13 '17 at 12:34

Be modern, use easiest modules for that (look at Problem intersecting shapefiles with OGR in Python)

from shapely.geometry import shape
import fiona
with fiona.open("ReBu.shp") as input:
    meta = input.meta
    with fiona.open('result.shp', 'w',**meta) as output:
       for feature in input:
         if abs(shape(feature['geometry']).area - 39.8953) < 0.00005:
             output.write(feature)

1) open the shapefile with Fiona: all the features are Python dictionaries
2) copy the properties (geometry, crs, fields definitions) of the original shapefile
3) create a new shapefile with these properties
4) transform the record geometry to a shapely geometry (for the area)
5) use the test of Vince
6) if True, write the record to the new shapefile

How to select polygons by criteria and save them to shp file?

2 Answers2