7

My objective is to create all possible pairs of connecting vertices in a concave polygon, i.e. only keeping those lines which are covered entirely by the polygon itself.

enter image description here

I have come up with the following query which is slow when it comes to more detailed and larger polygons. I am wondering if there are any handy tricks I could exploit to speed it up?

WITH poly AS (
    SELECT ST_GeomFromText('POLYGON((25.390624999999982 23.62298461759423,18.183593749999982 19.371888927008566,7.812499999999982 17.87273879517762,5.878906249999982 24.90497143578641,9.570312499999982 25.223427998254586,12.734374999999982 25.064303191014304,15.195312499999982 30.048744443788348,22.578124999999982 30.352588399664125,24.687499999999982 25.857838723065772,25.390624999999982 23.62298461759423))') AS geom
),
points AS (
    SELECT (g.gdump).path AS id, 
        (g.gdump).geom AS geom
    FROM (
        SELECT ST_DumpPoints(geom) AS gdump 
        FROM poly 
    ) AS g
)
SELECT 
ST_MakeLine(a.geom, b.geom) AS geom
FROM points a
CROSS JOIN points b
JOIN poly
ON ST_Contains(poly.geom, ST_MakeLine(a.geom, b.geom))
WHERE a.id < b.id;

From analyzing the query I guess that the JOIN ON ST_Contains is poorly implemented?

   Join Filter: st_contains(poly.geom, st_makeline(a.geom, b.geom))
   Rows Removed by Join Filter: 19
   CTE poly
     ->  Result  (cost=0.00..0.01 rows=1 width=32) (actual time=190.243..190.244 rows=1 loops=1)
   CTE points
     ->  Subquery Scan on g  (cost=0.00..40.03 rows=1000 width=64) (actual time=0.037..0.058 rows=10 loops=1)
           ->  ProjectSet  (cost=0.00..30.03 rows=1000 width=32) (actual time=0.034..0.049 rows=10 loops=1)
                 ->  CTE Scan on poly poly_1  (cost=0.00..0.02 rows=1 width=32) (actual time=0.000..0.002 rows=1 loops=1)
   ->  CTE Scan on poly  (cost=0.00..0.02 rows=1 width=32) (actual time=190.251..190.252 rows=1 loops=1)
   ->  Nested Loop  (cost=0.00..22530.00 rows=333333 width=64) (actual time=0.048..0.135 rows=45 loops=1)
         Join Filter: (a.id < b.id)
         Rows Removed by Join Filter: 55
         ->  CTE Scan on points a  (cost=0.00..20.00 rows=1000 width=64) (actual time=0.039..0.041 rows=10 loops=1)
         ->  CTE Scan on points b  (cost=0.00..20.00 rows=1000 width=64) (actual time=0.000..0.004 rows=10 loops=10)
 Planning Time: 0.290 ms
 JIT:
   Functions: 13
   Options: Inlining true, Optimization true, Expressions true, Deforming true
   Timing: Generation 3.830 ms, Inlining 18.323 ms, Optimization 121.388 ms, Emission 50.087 ms, Total 193.629 ms
 Execution Time: 196.074 ms
(21 rows)

Clarification:

There exist 2 issues at hand.

  1. For polygons without holes, a set of experiments were conducted with the suggestions made in the comments. I didn't notice any striking speed up.

  2. Polygons with holes have a negative effect on performance impact when using ST_Contains. Hence, the example provided is misleading for this investigation. Using the logic suggested in the accepted answer with the following example polygon consisting of approx. 2700 vertices yields a speed up of ~80% (20 minutes vs. 3 minutes). Polygon with hole

12
  • 3
    Why do you think its "poorly implemented"? Its the only part of that query that isn't simply moving numbers around, so it is bound to take the most time.
    – Spacedman
    Aug 10 at 10:57
  • 1
    Within/Contains performance is dependent on the efficiency of the index and the number of vertices in the polygon. Since the geometry is singleton, the index is moot. From there it's just N^2/2 pairs to test.
    – Vince
    Aug 10 at 11:59
  • 1
    I've tried to think of ways that might speed it up, for example by dividing the poly into convex subsets, or triangulating it. One algorithm might be to first work out what part of the polygon can be "seen" from the first vertex, and that might eliminate a lot of possible tests of things hidden round corners. I'm not sure how easy it would be to do in SQL though, and the resulting algorithm would still be O(n^2) for n vertices, and would be slower for convex polygons or not-very-concave ones.
    – Spacedman
    Aug 10 at 13:27
  • 1
    what about saving individual edge segments in an indexed column, then checking that your node to node line intersects exactly 4 edge segments?
    – JGH
    Aug 11 at 12:05
  • 1
    ah yes, they would have to be filtered out, like checking if the mid point (any point but the 2 ends) intersects with the source polygon - of why not a subdivided version of it.
    – JGH
    Aug 11 at 14:41
2

So, if I understand the question correctly, try to help the tools solve the problem you set, with ST_Contains().

What, you need to pay attention to:

  1. coordinate values have different lengths;
  2. polygons can have holes.

I tested the construction on the third polygon, which describes the coast "Chesapeake Bay", near Washington city (ST_GeomFromGeoJSON('{"type": "Polygon", "coordinates":[[[-75.7853058,37.100285],...), and in the above example left the example from your question, it is impossible to insert it into answer body :-).

Run the construction and check the result,

WITH poly AS 
(SELECT ST_MakePolygon(ST_ExteriorRing((ST_Dump(geom)).geom)) geom FROM
    (SELECT ST_GeomFromText('POLYGON((25.390624999999982 23.62298461759423,18.183593749999982 19.371888927008566,7.812499999999982 17.87273879517762,5.878906249999982 24.90497143578641,9.570312499999982 25.223427998254586,12.734374999999982 25.064303191014304,15.195312499999982 30.048744443788348,22.578124999999982 30.352588399664125,24.687499999999982 25.857838723065772,25.390624999999982 23.62298461759423))') AS geom
) boo),
poly_ext AS (SELECT ST_ExteriorRing((ST_Dump(geom)).geom) geom FROM poly),
points AS (SELECT (g.gdump).path AS id, (g.gdump).geom AS geom
    FROM (SELECT ST_DumpPoints(geom) AS gdump FROM poly_ext) AS g)
SELECT ST_MakeLine(a.geom, b.geom) AS geom FROM points a CROSS JOIN points b
JOIN poly ON ST_Contains(poly.geom, ST_MakeLine(a.geom, b.geom)) WHERE a.id < b.id;

Original spatial solutions...

Translated with www.DeepL.com/Translator (free version)

5
  • 1
    Thanks Cyril, I tested this and it is indeed faster - so basically the performance increase is coming from ST_ExteriorRing ? Can you elaborate a little more? Aug 11 at 20:16
  • PostGIS developers need to pay attention to the question you formulated, as I suggested a workaround to solve the problem, but the question itself needs careful consideration... Aug 14 at 20:36
  • 1
    I have added some additional context, I hope it makes more sense now. Aug 15 at 13:33
  • 1
    I think it doesn't hurt to be more specific, but the problems may be different... Aug 15 at 13:37
  • 1
    Agreed, I added more context now, thank you. Aug 15 at 13:43

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.