I have a shapefile that looks like this: enter image description here

the yellow polygone is a sqare of full extent size (not separate polygons) with class 0. Polygons of other colors have another classes (1, 2). If to look closer to the left side one can see small holes:

enter image description here

The same is in the upper left corner: enter image description here

Additionally, there are situations where 2 Polygons are touching each other and both touches the hole:

enter image description here

I need to join each hole that touches the boarder to the Polygon of classes 1 and 2 that touches that hole and boarder. At the same time I should not touch the "main part" of the class 0 (yellow) Polygon (I tried to highlight it by red color): enter image description here

The problem is that I can not come up with the feature for my python-code saing which "hole" of class 0 should belong to which Polygon of class 1 or 2. Seems like it should sound like "if the hole is surrounded by Polygon of class 1 or 2, join them", but I guess there is more determinant explaination.

  • What class will you assign to the polygons that share borders with class 1 and class 2?
    – Mayo
    May 22 at 14:02
  • it doesn't matter. May 23 at 6:06

1 Answer 1


I think the rule for joining extension holes with features of the main layer is: for every difference between the extension and the main layer (holes), which touches the extension and only touches 1 or two different classes of the main layer, must be joined with surrounding polygon from the main layer.

The way I've written this script makes it unnecessary to have the extension as a feature. So if you have the extension as a feature in your layer you need to remove it before running this script.

The implementation in GDAL looks something like this:

First import GDAL and set input and output layer.

from osgeo import gdalconst, ogr

input_layer_path = '' # set here the path of your input layer
input_driver_name = 'ESRI Shapefile' # set the file type of your input file
output_layer_path = '' # set here the path of your output layer
output_driver_name = 'ESRI Shapefile' # set the file type of your output file
class_field_name = '' # set here the name of the field that contains the classes
extent_class = 0 # set value to assign to holes in the middle of the layer
# get drivers for input and output layers
input_driver = ogr.GetDriverByName(input_driver_name)
output_driver = ogr.GetDriverByName(output_driver_name)
# get data sources for input and output layers
input_datasource = input_driver.Open(input_layer_path, gdalconst.GA_ReadOnly)
output_datasource = output_driver.CreateDataSource(output_layer_path)
# verify if data sources are valid
if input_datasource is None:
    raise FileNotFoundError('Input file not found')
if output_datasource is None:
    raise FileNotFoundError('Output file not found')
# get layers for input and output layers
input_layer = input_datasource.GetLayer()
output_layer = output_datasource.CopyLayer(input_layer, 'Output')

Then, get the difference between input layer and its extent.

# get the extent of the input layer
extent_points = input_layer.GetExtent()
# Create extent as a linear geometry
extent_line = ogr.Geometry(ogr.wkbLineString)
extent_line.AddPoint(extent_points[0], extent_points[3])
extent_line.AddPoint(extent_points[1], extent_points[3])
extent_line.AddPoint(extent_points[1], extent_points[2])
extent_line.AddPoint(extent_points[0], extent_points[2])
extent_line.AddPoint(extent_points[0], extent_points[3])
# Create polygon to compute the differences with input layer
extent_polygon = ogr.ForceToPolygon(extent_line)
# Get bounding line to filter only holes that are in the boundaries
bounding_line = extent_line
# Get difference between input layer and extent
holes = extent_polygon
for feature in input_layer:
    feature_geometry = feature.GetGeometryRef()
    holes = holes.Difference(feature_geometry)
hole_list = list(holes)

Now, hole_list is a list of polygons representing holes of the main layer. So now, we can filter all the holes to get only the ones that touch the extent and are touching only one class or two different classes.

# define function to get hole neighboring polygons from the input layer
def getNeighbors(geometry, layer):
    for feature in layer:
        feature_geometry = feature.GetGeometryRef()
        if feature_geometry.Touches(geometry):
            yield feature

feature_defn = output_layer.GetLayerDefn()
# loop through features of the main layer and join the holes to it
for feature in list(input_layer):
    feature_geometry = feature.GetGeometryRef()
    neighbor_holes = []
    # get a list of holes that are touching the feature
    for geometry in hole_list:
        if geometry.Touches(feature_geometry):
    output_geometry = ogr.Geometry(ogr.wkbPolygon)
    output_geometry = output_geometry.Union(feature_geometry)
    for hole in neighbor_holes:
        hole_neighbors = list(getNeighbors(hole, input_layer))
        # join if the hole is in the border and touches only one class or two different classes
        if hole.Touches(extent_line) and (len(hole_neighbors) == 1 or (len(hole_neighbors) == 2 and hole_neighbors[0][class_field_name] != hole_neighbors[1][class_field_name])):
            output_geometry = output_geometry.Union(hole)
        # else create a new different feature with the hole geometry
            new_feature = ogr.Feature(feature_defn)
            new_feature.SetField(class_field_name, extent_class)
        # delete the hole that was joined to avoid repeating the same geometry in next iterations
        del hole_list[hole_list.index(hole)]

And finally, save the output layer.

flush_err = output_datasource.SyncToDisk()
if not flush_err:
    print('Layer saved successfully')
    print('Error ' + str(flush_err))

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.