How to align edges of overlapping polygons in the middle line?

Question

The problem

Here are two polygons that overlap.

My ultimate goal is to reshape the grassland (green poly) and the reservoir (blue poly) so they share the edge somewhere in the middle of their intersection. It is acceptable if there would be a certain gap between the borders as long as they are not overlapping.

So, the area marked as W on the figure would be merged with the water body (blue poly to the left) and the area marked as G would be merged with the grassland (green poly to the right).

I have drawn the red line manually to illustrate the imaginary central line that splits the shared area into two equal (or roughly equal) pieces. The end result would look like this then (after merging):

Each polygon's shape was increased and they "are meeting in the middle".

I need to be able to do this analysis on any pair of overlapping polygons in an automated fashion. I am open to any existing solution. I am able to pull my data from an Esri geodatabase into other formats and do the processing somewhere else. If there are hints on the code libraries or tools for this kind of problem, I would be happy to learn more about them.

Research

There are some tools related to edge matching that exist in ArcGIS, mostly topology related. None of them supports splitting the overlapping polygons; they only allow choosing to what edge a polygon should be aligned to.

The Integrate tool would solve my problem in case the polygons are overlapping just over a minor area, but first the polygons can overlap over a significant area, and second the Integrate tool also modifies other vertices of polygons (not participating in the overlap) which I want to leave as is.

I got some crazy ideas about creating the longest central line over each polygon that would result when intersecting the water body and grassland, but I cannot really think through this using limited computational geometry knowledge I have. I guess the longest line that would split the polygon into two roughly equal parts would not always be the correct line I need to use.

I have started considering converting the polys to raster and reclassifying the overlapping pixels using an interpolation method (how far each pixel is from water and grass pixels), but can't figure out a workflow for doing this using ArcGIS Spatial Analyst tools.

Are there any good ideas on how to solve the problem?

Answer 1

1. Union the polygons
2. Convert to lines
3. Remove lines where FID_Reservoir OR FID_Grassland = -1 (no overlap)
4. Run 'Collapse Dual Lines to Centerline' tool
5. Use this tool (Split Polygons with Lines) to split the union result
6. Merge the polygons based on a common attribute (Reservoir ID; Grassland ID) using the dissolve tool (if you don't have a common attribute, add a unique id to the source data)

If you follow these steps, you should end up with the polygons split at the center of the overlap, as illustrated by your third image in your question.

Answer 2

Here is a code that uses centreline lable (https://github.com/ungarj/label_centerlines) library which can potentially solve your problem. Note this works pretty well with two polygons, but anything over three polygons becomes more difficult. The list of polygons consists of polygons in shapely Polygon format.

import geopandas as gdp 
from shapely.geometry import Polygon, MultiPolygon, MultiLineString
import fiona
from label_centerlines import get_centerline


for t in range(len(list_of_polygons) - 1, 0, -1):
# start loop through the other side of the polygon list
    for x in range(0, len(list_of_polygons) - 1):
        if list_of_polygons[x] != list_of_polygons[t]:
            if list_of_polygons[t].is_valid & list_of_polygons[x].is_valid:
                # boolean check to see if the polygons intersect
                if list_of_polygons[t].intersects(list_of_polygons[x]):
                    poly_list = []
                    intersection_list = []
                    intersection = list_of_polygons[t].intersection(list_of_polygons[x])
                    intersection_list.append(intersection)
                    if len(intersection_list[0]) > 1:
                        for i in intersection_list[0]:
                            if type(i) == type(list_of_polygons[0]):
                                poly_list.append(i)
                    else:
                        if type(intersection_list[0]) == type(list_of_polygons[0]):
                            poly_list.append(intersection_list[0])
                    for i in range(0, len(poly_list)):
                        # create centreline of the intersection
                        if (len(poly_list[i].exterior.coords[:]) > 20) & (poly_list[i].is_valid):
                            centreline = get_centerline(poly_list[i], segmentize_maxlen=0.000001, 
                                    max_points=6000, 
                                    simplification=0.000001,
                                    smooth_sigma=0.1) 
                            # centerline_list.append(centreline)
                            # change the geometry of current polygon
                            common = []
                            first = []
                            for y in range(0, len(list_of_polygons[x].exterior.coords[:])):
                                first.append(list_of_polygons[x].exterior.coords[y])
                                if list_of_polygons[x].exterior.coords[y] in poly_list[i].exterior.coords[:]:
                                    common.append(y)
                                else:
                                    continue
                            first[common[0]:common[-1]] = centreline.coords[:]
                            list_of_polygons[x] = Polygon(first)


                            # change the geometry of previous polygon
                            common = []
                            first = []
                            for y in range(0, len(list_of_polygons[t].exterior.coords[:])):
                                first.append(list_of_polygons[t].exterior.coords[y])
                                if list_of_polygons[t].exterior.coords[y] in poly_list[i].exterior.coords[:]:
                                    common.append(y)
                                else:
                                    continue
                            first[common[-1]:common[0]] = centreline.coords[:]
                            list_of_polygons[t] = Polygon(first)