9

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:

Non-fuzzy version

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.

  • I don't understand your operation. It seems you are looking for a common coarsening of two partitions. Without additional criteria, though, there will usually be many solutions. For instance, since coarsening (rather than refinement) appears to be your goal, why stop where you did? Why not just draw the common bounding square? – whuber Jul 23 '15 at 14:51
  • 1
    Thanks, I've mislabelled this. What I think I mean is finest common partition, or perhaps "least coarse". – EconAndrew Jul 23 '15 at 20:43
  • In that case the result would look very different from what you have drawn. It would be 4 x 4 chessboard of squares. From this one example I have been unable to deduce the rule you want to follow. Maybe you're trying to keep all edges common to all input features? What is the actual problem you are trying to solve? Could you provide a concrete example to help us understand what your question ought to be? – whuber Jul 23 '15 at 21:02
  • I've elaborated a lot - maybe this will help. It is true that in the fuzzy case I can't precisely specify my question, but I think in the exact case I know precisely what I mean (even if I'm not expressing it well). – EconAndrew Jul 23 '15 at 21:56
  • Thank you for those efforts (+1). In terms of your set-theoretic notation, partitions of a region form a partially ordered set: partition A is a refinement of B, and B is a coarsening of A, when every set in A is a subset of one in B. Your operation of combining appears to be the finest common coarsening of A and B. One way to tackle your fuzzy version would be to exploit the GIS's capabilities to remove dangles and slivers to correct small discrepancies between the two layers and then perform the non-fuzzy operation. – whuber Jul 23 '15 at 22:03
2

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

a b

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

symdiff

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

diff_symdiff

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.

  • Could you be explicit about how this algorithm either handles the fuzzy version of the problem being asked about, or else how it could be adapted to that version? – whuber Jul 26 '15 at 16:15
0

Errors free algorithm.

First set: enter image description here Second set: enter image description here

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:

enter image description here

Selected rows make your answer:

enter image description here

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

  • This looks like it will work correctly only in very special circumstances. How do you guarantee that you will end up with a non-overlapping, exhaustive partition of the common region? – whuber Jul 23 '15 at 14:52
  • Correct. Additional steps a)union datasets in terms of arcgis Union tool b) take first largest from merged table and check fraction of others inside c) remove others with fraction greater threshold, e.g. 90%. How is this? – FelixIP Jul 23 '15 at 18:59
  • I don't know, because I haven't figured out what the question is really asking yet. – whuber Jul 23 '15 at 19:00
  • Make up area using largest possible blocks. This is my understanding of question – FelixIP Jul 23 '15 at 19:02
  • The solution to that is to use a single block (the union of them all)! – whuber Jul 23 '15 at 19:17

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.