Removing polygons based on size of polygons in Python without using ArcPy?

Question

I am new to GIS. I want to remove some smaller polygons and only leave the polygons with a larger areas in Python scripts. Like this:

Here is the source files Shape files: buildings. I am thinking whether I can set a threshold for the size of the areas then remove them, but I don't know how to do this in Python.

Using the help by @Iridium , then I am able to compute the area of the polygons by this link Calculate area of polygons using OGR in python script

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile_to_be_used 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("Area", ogr.OFTReal)
new_field.SetWidth(32)
new_field.SetPrecision(2) #added line to set precision
layer.CreateField(new_field)

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

dataSource = None 

But how to filter out the polygons with small areas?

Answer 1

If you want a way to do these simple types of spatial operations in Python you can try geopandas which is a spatial extension to the (super popular and useful) pandas library.

As a quick example, firstly read the data in as a Geodataframe object, and reproject to a local projection (not being from South Africa, I hope this is right!).

import geopandas as gpd

buildings = gpd.read_file('buildings.shp').to_crs(epsg=22281)

At this point it's pretty easy to calculate the area for each geometry, and then threshold by some arbitrary number.

mask = buildings.area > 1000.  # metres squared

Finally, use the mask and fancy pandas indexing to mask out the points of interest:

selected_buildings = buildings.loc[mask]

At this point you can plot the data out:

ax = buildings.plot(color="gray", figsize=(20, 8))
selected_buildings.plot(ax=ax)

Or save the data back out to disk:

selected_buildings.to_file('filtered_buildings.shp')

Answer 2

Try using OGR/GDAL to achieve your task. One easy way to do this is calculate first the geometry in ArcGIS/QGIS (if possible), then iterate over the features with OGR and remove the ones with a value lower than X. Of course you could also calculate the area with GDAL/OGR like explained in this question.

A good place to start is the Python GDAL/OGR cookbook here

Answer 3

I hope it's not too late. I managed to do the same thing this way:

for feature in layer:
        geom = feature.GetGeometryRef()
        area = geom.GetArea() 
        print area
        if area<=320000:(or any area you want instead of 320000)
            layer.DeleteFeature(feature.GetFID())
            continue
        feature.SetField("Area", area)
        layer.SetFeature(feature)

I found a similar script here: https://www.programcreek.com/python/example/91667/osgeo.gdal.Polygonize

Hope it helps.