Python – How to Remove Small Holes from a Polygon Using GDAL and Shapely

gdalgeometrypolygonpythonshapely

I am writing a post-processing flow on Python. I am trying to come up with an algorithm to remove the small holes from the polygon.

I am using osgeo OGR right now, but I'm ok to do this part with shapely instead.

Currently, my solution is using Buffer like: geometry.Buffer(0.2 * 1 / 111111)

This works most of the time, but sometimes GDAL eats the whole memory on the machine and just crashes. I guess this is not optional at all. How can I accomplish my desired result?

The long gap on the right must persist, only the holes that are smaller the certain threshold have to be removed.

Best Answer

You can use the following example script. The script removes all holes from a polygon whose area is smaller than eps.

from shapely import wkt
from shapely.geometry import Polygon

# sample polygon
polygon_wkt = 'Polygon ((670651.58235156117007136 4087356.72129086637869477, 670941.20806862262543291 4087351.57990169199183583, 670944.42234629089944065 4087244.82082350691780448, 670648.40824284462723881 4087249.48037469061091542, 670651.58235156117007136 4087356.72129086637869477),(670669.30897117499262094 4087340.80642584012821317, 670671.69532915309537202 4087328.12050161883234978, 670709.06604568683542311 4087324.62686039274558425, 670709.12268763419706374 4087339.48419455718249083, 670669.30897117499262094 4087340.80642584012821317),(670670.63015638350043446 4087310.41313102655112743, 670670.24676115158945322 4087294.13430501008406281, 670711.11178680183365941 4087293.54062445322051644, 670711.46667919436004013 4087311.23376588197425008, 670670.63015638350043446 4087310.41313102655112743),(670669.74969205493107438 4087283.51269924966618419, 670668.85500870971009135 4087257.31934182625263929, 670755.23019460134673864 4087254.81062116287648678, 670754.28513388405553997 4087284.15052131982520223, 670669.74969205493107438 4087283.51269924966618419),(670796.18594406009651721 4087338.05096173053607345, 670797.80696410103701055 4087257.43525168858468533, 670807.4089543444570154 4087256.92088755592703819, 670806.12872222356963903 4087338.25089854281395674, 670796.18594406009651721 4087338.05096173053607345),(670821.05711140809580684 4087337.8436735225841403, 670824.12722124985884875 4087255.84220110531896353, 670929.67811466054990888 4087253.72090833634138107, 670930.47750282462220639 4087337.56890620524063706, 670821.05711140809580684 4087337.8436735225841403))'
polygon = wkt.loads(polygon_wkt)

list_interiors = []
eps = 1000

for interior in polygon.interiors:
    p = Polygon(interior)    
    if p.area > eps:
        list_interiors.append(interior)

new_polygon = Polygon(polygon.exterior.coords, holes=list_interiors)

Before:

After:

For a multipolygon, you should iterate over the list of its polygons as well.

from shapely import wkt
from shapely.geometry import Polygon, MultiPolygon

multipolygon_wkt = 'MultiPolygon (((670651.58235156117007136 4087356.72129086637869477, 670941.20806862262543291 4087351.57990169199183583, 670944.42234629089944065 4087244.82082350691780448, 670648.40824284462723881 4087249.48037469061091542, 670651.58235156117007136 4087356.72129086637869477),(670669.30897117499262094 4087340.80642584012821317, 670671.69532915309537202 4087328.12050161883234978, 670709.06604568683542311 4087324.62686039274558425, 670709.12268763419706374 4087339.48419455718249083, 670669.30897117499262094 4087340.80642584012821317),(670670.63015638350043446 4087310.41313102655112743, 670670.24676115158945322 4087294.13430501008406281, 670711.11178680183365941 4087293.54062445322051644, 670711.46667919436004013 4087311.23376588197425008, 670670.63015638350043446 4087310.41313102655112743),(670669.74969205493107438 4087283.51269924966618419, 670668.85500870971009135 4087257.31934182625263929, 670755.23019460134673864 4087254.81062116287648678, 670754.28513388405553997 4087284.15052131982520223, 670669.74969205493107438 4087283.51269924966618419),(670796.18594406009651721 4087338.05096173053607345, 670797.80696410103701055 4087257.43525168858468533, 670807.4089543444570154 4087256.92088755592703819, 670806.12872222356963903 4087338.25089854281395674, 670796.18594406009651721 4087338.05096173053607345),(670821.05711140809580684 4087337.8436735225841403, 670824.12722124985884875 4087255.84220110531896353, 670929.67811466054990888 4087253.72090833634138107, 670930.47750282462220639 4087337.56890620524063706, 670821.05711140809580684 4087337.8436735225841403)),((670648.88963292469270527 4087236.86040469771251082, 670943.25558728328906 4087232.40550666907802224, 670942.91290281957481056 4087129.60016754502430558, 670646.83352614217437804 4087133.71238110959529877, 670648.88963292469270527 4087236.86040469771251082),(670667.39459396700840443 4087217.67007472785189748, 670718.45457906532101333 4087218.35544365551322699, 670716.3984722828026861 4087206.36148742400109768, 670667.73727843072265387 4087207.04685635166242719, 670667.39459396700840443 4087217.67007472785189748),(670666.02385611203499138 4087192.65410887403413653, 670736.27417118020821363 4087192.65410887403413653, 670734.90343332523480058 4087158.38566249934956431, 670665.3384871844900772 4087159.07103142701089382, 670666.02385611203499138 4087192.65410887403413653),(670752.72302544000558555 4087176.54793907795101404, 670773.96946219238452613 4087175.17720122309401631, 670770.19993309106212109 4087158.72834696341305971, 670750.32423419377300888 4087157.70029357215389609, 670752.72302544000558555 4087176.54793907795101404),(670794.18784555338788778 4087175.86257015075534582, 670922.00915053102653474 4087176.20525461435317993, 670919.26767482107970864 4087139.19533252995461226, 670794.18784555338788778 4087140.90875484840944409, 670794.18784555338788778 4087175.86257015075534582),(670806.52448624826502055 4087215.61396794533357024, 670825.02944729058071971 4087214.92859901767224073, 670821.94528711691964418 4087184.08699728036299348, 670803.09764161077328026 4087183.40162835316732526, 670806.52448624826502055 4087215.61396794533357024),(670838.73682584054768085 4087214.24323009047657251, 670923.72257284971419722 4087213.55786116281524301, 670922.69451945857144892 4087182.37357496190816164, 670838.05145691300276667 4087186.48578852694481611, 670838.73682584054768085 4087214.24323009047657251)))'
multipolygon = wkt.loads(multipolygon_wkt)

list_parts = []
eps = 1000

for polygon in multipolygon.geoms:
    list_interiors = []
    
    for interior in polygon.interiors:
        p = Polygon(interior)

        if p.area > eps:
            list_interiors.append(interior)

    temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors)
    list_parts.append(temp_pol)
    
new_multipolygon = MultiPolygon(list_parts)

Before:

After:

Related Solutions

Remove Slivers in ArcGIS – Removing Small Spaces (Slivers) Between Polygons in ArcGIS

Use Eliminate (ArcInfo License required)

Eliminates polygons by merging them with neighboring polygons that have the largest area or the longest shared border. Eliminate is often used to remove small sliver polygons that are the result of overlay operations, such as Intersect or Union.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00170000005p000000.htm

enter image description here

ArcGIS 10 Users Note there is a NEW Feature in the Options for this Tool.

It is now possible to preserve the original POLYGON or POLYLINE with the option > ex_features (Optional) An input polyline or polygon feature class or layer that defines polygon boundaries, or portions thereof, that should not be eliminated.

Shapely – Removing Small Polygon Gaps in Polygons

The thing you are looking at is a sliver geometry. Similar to @sgillies's answer, except use a few buffer parameters to control the chiselled geometry shape:

import json
from shapely.geometry import shape, JOIN_STYLE
eps = 0.001  # epsilon that is approx. the width of slivers, e.g. 1 mm

# Load the original polygon from GeoJSON
poly = shape(json.loads('{"type": "Polygon", "coordinates": [[[...]]]}'))

# Here's the algorithm
fx = poly.buffer(eps, 1, join_style=JOIN_STYLE.mitre).buffer(-eps, 1, join_style=JOIN_STYLE.mitre)

# Compare number of vertices in the exterior LinearRing
print(len(poly.exterior.coords))  # 136
print(len(fx.exterior.coords))    # 135

Note that the fixed fx geometry has one less coordinate, which was the dangling sliver. Also note that some of the vertices may have wiggled around from their original position, usually several times less than eps.

Best Answer

Related Solutions

Remove Slivers in ArcGIS – Removing Small Spaces (Slivers) Between Polygons in ArcGIS

Shapely – Removing Small Polygon Gaps in Polygons

Related Question