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?
No comments:
Post a Comment