Finding truly intersected polygons

Question

Intro

Short version: My goal is to find geometries intersected by a target polygon.

Long version: We are building an interface for the user to upload new geometries (political maps) and allow him to solve 4D overlaps in a semi-automatic way. Steps: Find all intersecting polygons; Based on user review, modify polygons by using ST_Difference function.

To find overlaps - the most straightforward choice is to use ST_Overlaps or ST_Intersects. But this method will give us an unlimited amount of false-positive results. Probably because of the tolerance, borders located next to each other would be treated as an intersection.

SELECT *, ST_Intersects(
    ( SELECT territory FROM diff_test WHERE id = 1 ),
    ( territory )
) as intersects
FROM diff_test
WHERE diff_test.id != 1

The other method I can think of is to select every polygon limited by bounding box. But then there is a need to filter out all the "touching" borders, and ST_Touching is not helpful if you don't have common points. It will also not work.

Back to ST_Intersects:

Let's say we have two overlapped polygons - id 1 (1430) and id 2 (1487).

We've calculated two new polygons.

• id 3 (1430 - 1487) ST_Difference(id1, id2)
• id 4 (1487 - 1430) ST_Difference(id2, id1) # this one can be ignored

What did I expect:

• id 3 is inside id 1 (ST_Overlaps = false)
• id 3 is touching id 2 (ST_Touches = true)
• id 3 does not intersect with id 2 (ST_Intersects = false)

Nothing from this was as I've expected.

The query was executed targeting id 3

Problem and Questions

• What is the best solution to find truly intersect polygons?
• What would be the best practice of avoiding intersections from newly created polygon by ST_Difference and his counterpart?

Thoughts

• I can think of multiple ST_Difference operations, but it will not guarantee anything and it is slow. Next query still returns true instead of false

SELECT ST_Overlaps(t1mask, origt2) FROM ( SELECT
    (SELECT territory FROM diff_test WHERE id = 2) as origt2
    ,(SELECT ST_Difference(mask, intersects) FROM ( 
        SELECT
        ( SELECT territory FROM diff_test WHERE id = 3) as mask
        ,( 
            SELECT ST_Intersection(
            ( SELECT territory FROM diff_test WHERE id = 3 )
            ,( SELECT territory FROM diff_test WHERE id = 2 ))
        ) as intersects
    ) as t ) as t1mask
) as t

• As another option - cut a little bit more, to prevent this from happening.

• Or maybe this problem can be reproduced only on my data and can be solved by ST_SnapToGrid or something similar

Test sample

Link for this shapes

Answer 1

Uf you want to find overlaps of two polygons, then your best fried is st_intersection in conjunction with st_intersects (it's not the same function). You should firstly join polygons using st_intersects. Then get intersected geometries with st_intersection. This will produce different geometry collections based on type of interaction. Points and Lines in the case of touch and polygons in the case of true overlap. So you can then pretty simple pick what you want. BTW PostGIS performs lanient geometry marching, thus there shouldn't be requirements of common points to recognize spatial relation of two geometries...

If you want implement some sort of auto-complete polygon function, then I would do this:

1. Collect polygons from user
2. Find intersection of all existing polygons as described above
3. Union all intersected geometris to one geometry
4. Find difference between unioned and original geometry using st_difference.

I thing this should produce want you want...

Answer 2

To find polygons which "truly" intersect (ie. have some shared area in common) then ST_Relate must be used with an Intersection Matrix pattern of T********. There is no standard named predicate which provides this semantics. (It would be nice if there was a named predicate called something like "interiorIntersects".) The PostGIS doc has an example.

Another option is to compute the actual intersection geometry (using ST_Intersection) and check whether it has dimension = 2. However, this will be slower than using the predicate. Also, the intersection computation is subject to numerical round-off, which means that in certain cases the intersection may not contain an area even if the overlaps predicate returns true.

Note also that if the input data is not a true polygonal coverage, it may be that "adjacent" polygons actually have a non-zero intersection area (and hence ST_Relate(a, b, 'T********') = true). In this case it is necessary to compute the intersection area and then filter out cases with a small tolerance value on the intersection area.

Answer 3

I think it is likely that the analysis of ST_Overlaps and ST_Touches for multi-part geometries may not be entirely correct in versions of PostGIS less than 3.0.0. Not that of ST_Intersects, a function that supports geometry collections since version 2.5.0.

The product of the intersection between two geometries can result in a GEOMETRYCOLLECTION even if the original geometries were single-part polygons. This collection may contain polygons, lines and points, depending on the spatial condition between intersected geometries.

If what you want to obtain those intersections in which at least one of its parts is a polygon, you can filter the results so that they satisfy that ST_Dimension of the intersection returns the value 2:

SELECT *
FROM (
  SELECT a.id a_id, b.id b_id, ST_Intersection( a.territory, b.territory) inters_geom
  FROM t1 a
  INNER JOIN t1 b
    ON (a.id = 3
    AND ST_Intersects( a.territory, b.territory)))
WHERE ST_Dimension( inters_geom) = 2;  

---

A much more detailed analysis about the spatial condition of the analyzed geometries can be performed through ST_Relate (PostGIS >= 3.0.0 for GEOMETRYCOLLECTION support).

Answer 4

Thank you all for all the hints and advice.

So far I've ended up with filtering out polygons smaller than 100sqm considering them as an error. UPD: As I am storing geometry in non-projected way with SRID:4326 to calculate area of the polygon I am casting geometry to geography

SELECT geom, ST_Area(geom) sqft FROM (
    SELECT
        (ST_Dump(ST_Intersection(orig, diff))).geom as geom
    FROM (
        SELECT *
            , poly2 as orig, ST_Difference(poly1, poly2) as diff
--          , poly1 as orig, ST_Difference(poly2, poly1) as diff
        FROM (
            SELECT
                (SELECT territory FROM api_spacetimevolume WHERE id = 225) as poly1
                , (SELECT territory FROM api_spacetimevolume WHERE id = 230) as poly2
        ) as foo
    ) as foo
) as foo
WHERE ST_Dimension(geom) = 2 AND ST_Area(geom::geography) > 100

I am terrified in amount of time for this procedure.

• Find all candidates for intersection in space and time with ST_Intersects
• Filter out false-positive with function above
• Do it again ignoring some of the polygons
• Save the results

And it is not all, at some point there will be a necessity to 'heal' the data set from 'errors' introduced by ST_Difference

UPD, final version of the query

SELECT id, entity_id, end_date, start_date, ST_Union(xing) FROM (
    SELECT *,
      (ST_Dump(ST_Intersection(territory, diff))).geom as xing
    FROM (
      SELECT *,
        ST_Difference(
          territory,            
          'SRID=4326;POLYGON ((-132.890625 72.50172235139388, -115.3125 -55.77657301866769, 125.15625 -53.5403073915002, 120.234375 74.01954331150228, -132.890625 72.50172235139388))'::geometry
        ) as diff
      FROM (
        SELECT *
        FROM api_spacetimevolume as stv
        WHERE stv.end_date >= '0'::numeric(10,1) AND stv.start_date <= '999999999'::numeric(10,1)
        AND ST_Intersects(
          territory,
          'SRID=4326;POLYGON ((-132.890625 72.50172235139388, -115.3125 -55.77657301866769, 125.15625 -53.5403073915002, 120.234375 74.01954331150228, -132.890625 72.50172235139388))'::geometry)
      ) as foo
    ) as foo
  ) as foo
  WHERE ST_Dimension(xing) = 2 AND ST_Area(xing::geography) > 100
  GROUP BY id, entity_id, start_date, end_date