5

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

Example polygons

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.

Broken illusions 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

  • 1
    The only way to get adjacent polygons without overlaps or gaps is to make them use exactly same vertices on common boundaries which are curved. I suppose your false positives are rather true positives but overlaps are neglicible in practical use. Do you have a possibilty to fix the source data? – user30184 Nov 14 '19 at 9:11
  • I did the ST_Difference operation, which was intended to fix any overlaps, so it is truly false-positive. Even after ST_Difference ST_Intersects is returning true. Of course, I am applying ST_Intersects to a polygon created as a result of ST_Difference – MiklerGM Nov 14 '19 at 15:29
  • 1
    You could have a try by pushing the polygons from ST_Difference through ST_SnapToGrid when checking them wth ST_Intersects. It would also be interesting to get some test data from a couple of cases either as WKT or as a download. – user30184 Nov 14 '19 at 18:08
  • Sure, I've uploaded the diff_test table as is. I can publish more polygons tomorrow. ST_SnapToGrid looks like a good idea. I am thinking of processing the whole dataset with it even before ST_Difference – MiklerGM Nov 14 '19 at 22:33
  • 1
    All geometries intersects each other, that makes sense to me. If you are waiting a positive intersection when it returns a polygon, you can filter the intersections with dimension 2. Analyzing multi-part geometries will not be a problem? What version of PostGIS are you using? – Gabriel De Luca Nov 15 '19 at 6:59
3

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

  • It is not exactly the problem we are trying to solve. The problem is, even If I will apply ST_Difference(Poly1, Poly2) as Poly3. I will receive a non zero geometry in ST_Intersection(ST_Difference(Poly1, Poly2), Poly2) as intersection. And even after doing ST_Difference(Poly3, intersection), polygons would be marked as overlapped – MiklerGM Nov 14 '19 at 16:27
  • @MiklerGM st_intersects returns true If a geometry or geography shares any portion of space, even if they only touches. That's it. If you intersect touching geometries then you should get line representing common border, thus non null geometry. – DavidP Nov 14 '19 at 18:12
  • Well, I am trying to find a solution to avoid sharing territory between two polygons. I am using ST_Difference to subtract one from another, and it does not help. Btw, ST_Touches is false that means it is more than sharing the same border – MiklerGM Nov 14 '19 at 22:37
3

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.

  • Territory 3 is inside territory 1, so their intersection dimension is 2 but territory 3 doesn't overlap territory 1. In fact, the definition provided doesn't state that two polygons overlap if their intersection has dimension 2. Beyond that, if the geometry returned by the intersection between the two original geometries is not a GEOMETRYCOLLECTION, I agree that the error in PostGIS 2.5.0 may be due to the rounding of the coordinates of some vertices of the resulting geometry. – Gabriel De Luca Nov 17 '19 at 18:19
  • 1
    @GabrielDeLuca you're right, my mistake. ST_Overlaps doesn't detect wholly-contained cases. In fact, there is no named predicate which does exactly this. Instead, ST_Relate with a specify pattern of 'T********' (interior intersects) must be used. I'll edit my reply. – dr_jts Nov 17 '19 at 23:33
  • Territory 3 and 4 were produced in a result of ST_DIfference(1,2) and ST_Difference(2,1). And the problem is, that they intersects with parents. I guess that the problem lays somewhere in ST_Differnece, even with strong snap to grid ST_CollectionExtract( ST_Force_Collection( ST_MakeValid( ST_SnapToGrid(first_poly, 0.01)), 3)) – MiklerGM Nov 18 '19 at 18:45
  • 1
    Yes, in general the overlay operations won't produce results that "match" exactly with the inputs. The problem is that it's not possible to exactly represent a point which is the intersection of two line segments. in general. – dr_jts Nov 18 '19 at 21:02
0

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

  • The spatial predicates have been extensively tested in all versions of PostGIS, so they are correct. In this case it is how they are being used which is the issue. – dr_jts Nov 17 '19 at 17:43
0

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

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.