I am polygonizing a single band raster.


I am using gdal.Polygonize() function to vectorize the raster.

img= r"C:\Users\aoi2.tif"
src_ds = gdal.Open(img)
srcband = src_ds.GetRasterBand(1)

# create output datasource
dst_layername = "output"

drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername +".shp")
src = osr.SpatialReference(wkt=prj)

dst_layer = dst_ds.CreateLayer(dst_layername, srs = src)

gdal.Polygonize( srcband, srcband, dst_layer, -1, [], callback=None )

del src_ds, srcband, dst_ds, dst_layer

output shapefile looks like:


My aim is to erode the islands of pixels smaller than 3 pixels and keep the rest, how can I do that?


Two options for me :

  • You can do it directly on your raster using connectedComponentsWithStats of opencv and CC_STAT_AREA option. See this link (here).

  • You can do it on your shapefile. First, you have to calculate a new field "area" for your shapefile, then, you can do an attribute query and delete features under a specific value (area = 3 pixels in your case).

In this case, here, a function to calculate the polygon area in meters:

def polygon_area(shp):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    src = driver.Open(shp, 1)
    layer = src.GetLayer()

    new_field = ogr.FieldDefn("Area", ogr.OFTReal)

    for feature in layer:
        geom = feature.GetGeometryRef()
        area = geom.GetArea()
        feature.SetField("Area", area)

Here, a function to delete the features under a specific value:

def sort_by_area(shp, area):
    layer = shp.GetLayer()

    for feature in layer:
        shp.ExecuteSQL('REPACK ' + layer.GetName())
| 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.