Search code examples
postgresqlpostgis

Postgis : Query to get list of multiple neighbour polygon which area sum equal to 42

I'm finally asking a question for this because it's been days I'm on it!

My problem is the following one : In my database, I have a lot of polygons which are in the same region. I'm trying to create a query for 3 polygons which let me know all possibilities where 3 polygons are neighbour AND the sum of their 3 areas is superior to 41 and inferior to 43 (example)

I already got it for 2 polygons but all my tries for 3 polygons have been a failure, I can get a result but I have to wait way too long (a few hours sometimes ). However, with 2 polygons it can be a few seconds ... I'm running my script on my server (12 core and 64 RAM)

Here is the code I was using for 3 polygons search ...

SELECT DISTINCT
    p1.id As p1_id, p2.id As p2_id, p3.id As p3_id
FROM table1 As p1, table1 As p2, table1 As p3
WHERE 
    p1.region_id IN (17) 
    AND p2.region_id IN (17) 
    AND p3.region_id IN (17) 
    AND p1.id <> p2.id
    AND p1.id <> p3.id
    AND p2.id <> p3.id 
    AND ((ST_Area(p1.polygon) * 0.67 + ST_Area(p2.polygon) * 0.67 + ST_Area(p3.polygon) * 0.67) > 41.0)
    AND ((ST_Area(p1.polygon) * 0.67 + ST_Area(p2.polygon) * 0.67 + ST_Area(p3.polygon) * 0.67) < 43.0)
    AND (
    (ST_DWithin(p1.polygon, p2.polygon, 1.0, true) AND ST_DWithin(p2.polygon, p3.polygon, 1.0, true))
      OR (ST_DWithin(p1.polygon, p3.polygon, 1.0, true) AND ST_DWithin(p3.polygon, p2.polygon, 1.0, true)) 
      OR (ST_DWithin(p2.polygon, p1.polygon, 1.0, true) AND ST_DWithin(p1.polygon, p3.polygon, 1.0, true))
     )
LIMIT 10;

Any advice or help is welcome!

Another info: my goal was to be able to generate this kind of query for X number of polygons

Thank's and have a nice day

################################## Edit

As asked by Jim, here is an example

We suppose this region with 11 polygons:

We are looking for 3 polygons which are neighbour and sum area equal to 7 In this case, I suppose 3 solutions exist:

  1. P2-P3-P11
  2. P3-P11-p10
  3. P4-P10-P11

For the last one, it's only an intersect on the point, I'm not sure I will keep this result later, it will maybe 2 points at least to be considered (by this I'm avoiding angle ..)

enter image description here

P1 = POLYGON(( 0.0 0.0, 0.0 1.0, 1.0 1.0, 1.0 0.0 ))
P2 = POLYGON(( 1.0 0.0, 1.0 1.0, 3.0 1.0, 3.0 0.0 ))
P3 = POLYGON(( 3.0 0.0, 3.0 2.0, 4.0 2.0, 4.0 0.0 ))
P4 = POLYGON(( 0.0 1.0, 0.0 3.0, 1.0 3.0, 1.0 1.0 ))
P5 = POLYGON(( 1.0 1.0, 1.0 2.0, 2.0 2.0, 2.0 1.0 ))
P6 = POLYGON(( 2.0 1.0, 2.0 2.0, 3.0 2.0, 3.0 1.0 ))
P7 = POLYGON(( 1.0 2.0, 1.0 3.0, 2.0 3.0, 2.0 2.0 ))
P8 = POLYGON(( 2.0 2.0, 2.0 3.0, 3.0 3.0, 3.0 2.0 ))
P9 = POLYGON(( 0.0 3.0, 0.0 4.0, 1.0 4.0, 1.0 3.0 ))
P10 = POLYGON(( 1.0 3.0, 1.0 4.0, 3.0 4.0, 3.0 3.0 ))
P11 = POLYGON(( 3.0 2.0, 3.0 5.0, 4.0 5.0, 4.0 2.0 ))

So for a multipolygon:

MULTIPOLYGON((( 0.0 0.0, 0.0 1.0, 1.0 1.0, 1.0 0.0 )),
(( 1.0 0.0, 1.0 1.0, 3.0 1.0, 3.0 0.0 )),
(( 3.0 0.0, 3.0 2.0, 4.0 2.0, 4.0 0.0 )),
(( 0.0 1.0, 0.0 3.0, 1.0 3.0, 1.0 1.0 )),
(( 1.0 1.0, 1.0 2.0, 2.0 2.0, 2.0 1.0 )),
(( 2.0 1.0, 2.0 2.0, 3.0 2.0, 3.0 1.0 )),
(( 1.0 2.0, 1.0 3.0, 2.0 3.0, 2.0 2.0 )),
(( 2.0 2.0, 2.0 3.0, 3.0 3.0, 3.0 2.0 )),
(( 0.0 3.0, 0.0 4.0, 1.0 4.0, 1.0 3.0 )),
(( 1.0 3.0, 1.0 4.0, 3.0 4.0, 3.0 3.0 )),
(( 3.0 2.0, 3.0 5.0, 4.0 5.0, 4.0 2.0 )))

I hope I haven't made a mistake on coordinates. Thank's everyone for the differents answer :D

Solution

  • For 3, we can say that a central polygon is connected to 2 others. I would focus on finding this central polygon. It will eventually lead to duplicates (if the 3 polygons touch each others) that you would have to filter out.

    Also make sure to have an index on st_area(geom) or to pre-compute it. (In fact, a multi-column index on id, region, st_area(geom) should be the most efficient)

    Most importantly, have a spatial index.

    SELECT
    p1.id As p1_id, p2.id As p2_id, p3.id As p3_id
FROM table1 As p1
 JOIN table1 As p2 ON ST_DWithin(p1.polygon, p2.polygon, 1.0, true)
 JOIN table1 As p3 ON ST_DWithin(p1.polygon, p3.polygon, 1.0, true)
WHERE 
    p1.region_id IN (17) 
    AND p2.region_id IN (17) 
    AND p3.region_id IN (17) 
    AND p1.id <> p2.id
    AND p1.id <> p3.id
    AND p2.id <> p3.id 
    AND ((ST_Area(p1.polygon) * 0.67 + ST_Area(p2.polygon) * 0.67 + ST_Area(p3.polygon) * 0.67) > 41.0)
    AND ((ST_Area(p1.polygon) * 0.67 + ST_Area(p2.polygon) * 0.67 + ST_Area(p3.polygon) * 0.67) < 43.0
LIMIT 10;