Fuzzy least-coarse common partition algorithm

Question

Given two different partitions of a shape (for argument's sake, two different administrative divisions of a country), how can I find a new partition into which both of those partitions fit, allowing for (and optimising) some error?

For example, ignoring the error, I want an algorithm that does this:

Perhaps it helps to express this in set terms. Using the following numbering:

I can express the partitions above as:

A = {{1},{2},{3,4,7,8},{5},{6},{9,10,13,14},{11},{12},{15},{16}}

B = {{1,2,5,6},{3},{4},{7},{8},{9},{10},{13},{14},{11,15},{12,16}}

A dot B = {{1,2,5,6},{3,4,7,8},{9,10,13,14},{11,15},{12,16}}

and the algorithm for producing A dot B seems straightforward (something like, if two elements are in a partition together in A (B) merge the partitions they are in in B (A) - repeat until A and B are equal).

But now imagine that some of these lines are slightly different between the two partitions, so that this perfect answer is not possible, and instead I want the optimal answer subject to minimizing some error criterion.

Take a new example:

Here in the left column we have two partitions with no common lines (apart from the outside border itself). The only possible solution of the above kind is the trivial one, the right column. But if we allow "fuzzy" solutions, then the middle column may be permissible, with say 5% of the total area being contested (i.e. allocated to a different subarea in each coarsened partition). So we might describe the middle column as representing the "least-coarse common partition with <= 5% error".

Whether the actual answer is then the partition in top row, middle column or middle row, middle column - or something in between, is less important.

Answer 1

You can do this by evaluating the difference of a polygon's boundary to the symmetric difference between their boundaries, or symbolically expressed as:

Difference(a, SymDifference(a, b))

Take geometries a and b, expressed as MultiLinestrings over the next two lines and images:

MULTILINESTRING((0 300,50 300,50 250,0 250,0 300),(50 300,100 300,100 250,50 250,50 300),(0 250,50 250,50 200,0 200,0 250),(50 250,100 250,100 200,50 200,50 250),(100 300,200 300,200 200,100 200,100 300),(0 200,100 200,100 100,0 100,0 200),(100 200,150 200,150 150,100 150,100 200),(150 200,200 200,200 150,150 150,150 200),(100 150,150 150,150 100,100 100,100 150),(150 150,200 150,200 100,150 100,150 150))
MULTILINESTRING((0 300,100 300,100 200,0 200,0 300),(100 300,150 300,150 250,100 250,100 300),(150 300,200 300,200 250,150 250,150 300),(100 250,150 250,150 200,100 200,100 250),(150 250,200 250,200 200,150 200,150 250),(0 200,50 200,50 150,0 150,0 200),(50 200,100 200,100 150,50 150,50 200),(0 150,50 150,50 100,0 100,0 150),(50 150,100 150,100 100,50 100,50 150),(100 200,150 200,150 100,100 100,100 200),(150 200,200 200,200 100,150 100,150 200))

The symmetric difference, where portions of a and b do not intersect, is:

MULTILINESTRING((50 300,50 250),(50 250,0 250),(100 250,50 250),(50 250,50 200),(150 150,100 150),(200 150,150 150),(150 300,150 250),(150 250,100 250),(200 250,150 250),(150 250,150 200),(50 200,50 150),(50 150,0 150),(100 150,50 150),(50 150,50 100))

And finally, evaluate the difference between either a or b and the symmetric difference:

MULTILINESTRING((0 300,50 300),(0 250,0 300),(50 300,100 300),(100 300,100 250),(50 200,0 200),(0 200,0 250),(100 250,100 200),(100 200,50 200),(100 300,150 300),(150 300,200 300,200 250),(200 250,200 200),(200 200,150 200),(150 200,100 200),(100 200,100 150),(100 150,100 100),(100 100,50 100),(50 100,0 100,0 150),(0 150,0 200),(150 200,150 150),(200 200,200 150),(150 150,150 100),(150 100,100 100),(200 150,200 100,150 100))

You can implement this logic in GEOS (Shapely, PostGIS, etc.), JTS, and others. Note that if the input geometries are polygons, then their boundaries need to be extracted, and the result can be polygonized. For example, shown with PostGIS, take two MultiPolygons, and get a MultiPolygon result:

SELECT
  ST_AsText(ST_CollectionHomogenize(ST_Polygonize(
    ST_Difference(ST_Boundary(A), ST_SymDifference(ST_Boundary(A), ST_Boundary(B)))
  ))) AS result
FROM (
  SELECT 'MULTIPOLYGON(((0 300,50 300,50 250,0 250,0 300)),((50 300,100 300,100 250,50 250,50 300)),((0 250,50 250,50 200,0 200,0 250)),((50 250,100 250,100 200,50 200,50 250)),((100 300,200 300,200 200,100 200,100 300)),((0 200,100 200,100 100,0 100,0 200)),((100 200,150 200,150 150,100 150,100 200)),((150 200,200 200,200 150,150 150,150 200)),((100 150,150 150,150 100,100 100,100 150)),((150 150,200 150,200 100,150 100,150 150)))'::geometry AS a,
    'MULTIPOLYGON(((0 300,100 300,100 200,0 200,0 300)),((100 300,150 300,150 250,100 250,100 300)),((150 300,200 300,200 250,150 250,150 300)),((100 250,150 250,150 200,100 200,100 250)),((150 250,200 250,200 200,150 200,150 250)),((0 200,50 200,50 150,0 150,0 200)),((50 200,100 200,100 150,50 150,50 200)),((0 150,50 150,50 100,0 100,0 150)),((50 150,100 150,100 100,50 100,50 150)),((100 200,150 200,150 100,100 100,100 200)),((150 200,200 200,200 100,150 100,150 200)))'::geometry AS b
) AS f;
                               result
--------------------------------------------------------------------------------
MULTIPOLYGON(((0 300,50 300,100 300,100 250,100 200,50 200,0 200,0 250,0 300)),((100 250,100 300,150 300,200 300,200 250,200 200,150 200,100 200,100 250)),((0 200,50 200,100 200,100 150,100 100,50 100,0 100,0 150,0 200)),((150 200,200 200,200 150,200 100,150 100,150 150,150 200)),((100 200,150 200,150 150,150 100,100 100,100 150,100 200)))

Note that I have not extensively tested this method, so take these as ideas as a starting point.

Answer 2

Errors free algorithm.

First set: Second set:

Merge 2 sets and sort in descending order by area. Select rows in the table (top => down) until total of areas = total area (16 in this case) reached:

Selected rows make your answer:

The criteria is going to be a difference between accumulated areas and actual total.