0

I have polygonized a single band raster and created an area field. I want to keep the 10 largest polygons by area and delete the rest. The following script computes the area for all polygons and sorts the area field and keeps 10 largest values and assign NULL to rest of the features.

img = r"C:\Users\test\input_img.tif"   

# mapping between gdal type and ogr field type
type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                gdal.GDT_UInt16: ogr.OFTInteger,
                gdal.GDT_Int16: ogr.OFTInteger,
                gdal.GDT_UInt32: ogr.OFTInteger,
                gdal.GDT_Int32: ogr.OFTInteger,
                gdal.GDT_Float32: ogr.OFTReal,
                gdal.GDT_Float64: ogr.OFTReal,
                gdal.GDT_CInt16: ogr.OFTInteger,
                gdal.GDT_CInt32: ogr.OFTInteger,
                gdal.GDT_CFloat32: ogr.OFTReal,
                gdal.GDT_CFloat64: ogr.OFTReal}

ds = gdal.Open(img)
prj = ds.GetProjection()
srcband = ds.GetRasterBand(1)

# Create shapefile layer
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(r"C:\Users\test")
srs = osr.SpatialReference(wkt=prj)
dst_layer = dst_ds.CreateLayer("output_shp", srs=srs)

# Create area field
col = ogr.FieldDefn('AREA', ogr.OFTReal)
dst_layer.CreateField(col)

# Compute area and sort values    
features = []

def get_area(feature):
        area= feature.GetGeometryRef().GetArea()
        return area


for feat in dst_layer:
    features.append(feat)

# Sort 10 largest values
sorted_areaVals = sorted(features, key= get_area, reverse=True)[:10]

# Write sorted values to the filed
for sv in sorted_areaVals:

    sv.SetField("AREA", get_area(sv))
    dst_layer.SetFeature(sv)

Area field gets the 10 largest values while rest of features get NULL. How can I delete the features with NULL in area field?

1

If you want to keep the 10 largest polygons, I'd recommend executing a SQL sentence on your original data source and then copying the result to a new data source.

import ogr

ds = ogr.Open('some_file.shp')
lyr = ds.GetLayer()
name = lyr.GetName()
n = 10  # features to keep
sql = f'SELECT * FROM {name} ORDER BY OGR_GEOM_AREA DESC LIMIT {n}'  # this query selects the n largest features

lyr = ds.ExecuteSQL(sql)  # here you are storing the results of the query

# create new shapefile to store
driver = ogr.GetDriverByName('Esri Shapefile')
new_ds = driver.CreateDataSource('some_other_file.shp')
new_ds.CopyLayer(lyr, name)

del ds, lyr, new_ds

Note that OGR_GEOM_AREA is a special field that returns the feature's geometry area. You can read more information here. Also, this script will only work on Python versions above (and including) 3.6 as it uses f-strings. If you are using a lower Python version change the sql query line to:

sql = 'SELECT * FROM {} ORDER BY OGR_GEOM_AREA DESC LIMIT {}'.format(name, n)
| improve this answer | |

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.