Holes closing near extent boarder

Question

I have a shapefile that looks like this:

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:

The same is in the upper left corner:

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

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):

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.

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):
            neighbor_holes.append(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
        else:
            new_feature = ogr.Feature(feature_defn)
            new_feature.SetField(class_field_name, extent_class)
            new_feature.SetGeometry(hole)
            output_layer.CreateFeature(new_feature)
        # delete the hole that was joined to avoid repeating the same geometry in next iterations
        del hole_list[hole_list.index(hole)]
    feature.SetGeometry(output_geometry)
    output_layer.CreateFeature(feature)

And finally, save the output layer.

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