5

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:

enter image description here

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?

7

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)

Subsetted buildings in blue

Or save the data back out to disk:

selected_buildings.to_file('filtered_buildings.shp')
| improve this answer | |
1

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

| improve this answer | |
  • I have to use python to realize this and no Arcpy. I wish I could have some python scripts example. – statistics_learning Apr 27 '18 at 9:48
  • I edited my answer adding a link to calculate area of polygons with python and OGR. Shouldn't be complicated to edit the code to fit your needs. The Python GDAL/OGR cookbook I also linked you has tons of python scripts examples that directly fit your needs. – Iridium Apr 27 '18 at 9:53
  • Thank you so much. It is very helpful. But I still can't not find the solution to filter out small areas based on the size of the areas. I will try to do this. But if you have solution to this, I will appreciate very much. – statistics_learning Apr 27 '18 at 10:11
1

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.

| 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.