Postgis：查询以获取面积总和等于42的多个相邻多边形的列表

Question

我终于要问这个问题了，因为我已经好几天了！

我的问题如下：在我的数据库中，我有很多位于同一区域的多边形。我正在尝试为 3 个多边形创建一个查询，它让我知道 3 个多边形相邻并且它们的 3 个面积之和优于 41 且低于 43 的所有可能性（示例）

我已经获得了 2 个多边形的结果，但是我对 3 个多边形的所有尝试都失败了，我可以得到结果，但我必须等待太久（有时需要几个小时）。但是，对于 2 个多边形，可能需要几秒钟...... 我在我的服务器上运行我的脚本（12 核和 64 内存）

这是我用于 3 个多边形搜索的代码...

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;

欢迎任何建议或帮助！

另一个信息：我的目标是能够为 X 个多边形生成这种查询

谢谢，祝你有愉快的一天

################################## 编辑

正如吉姆所问，这里有一个例子

我们假设这个区域有 11 个多边形：

我们正在寻找 3 个相邻且总面积等于 7 的多边形在这种情况下，我假设存在 3 个解决方案：

P2-P3-P11
P3-P11-p10
P4-P10-P11

最后一个，只是一个点的交集，我不确定以后会不会保留这个结果，可能至少要考虑2个点（这样我就避免了角度.. )

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

所以对于多面体：

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

希望我没有弄错坐标。谢谢大家的不同回答 :D

Answer 1

对于 3，我们可以说一个中心多边形与另外 2 个多边形相连。我会专注于找到这个中心多边形。它最终会导致您必须过滤掉的重复项（如果 3 个多边形相互接触）。

还要确保在 st_area(geom) 上有一个索引或预先计算它。（其实id, region, st_area(geom)上的多列索引应该效率最高）

最重要的是，有一个空间索引。

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;

Answer 2

我想我有一个可扩展的解决方案...如果我们合并接触的多边形，我们只会遇到比较两个多边形的情况。如果我们知道哪些多边形接触，我们就可以在多边形 ID 数组中跟踪接触的多边形。

我使用 intarray 对 id 进行排序：

Create extension intarray;
alter table polygons add id_array  int[];
update polygons set id_array = array[id];
create index polygons_geom on polygons using gist(geom);

编辑 - 使用 select distinct on 而不是 group by，并基于 @JGH 的观点，我们可以进行一次空间计算，然后在数组中查找

预先计算出接触多边形和面积，就不需要再进行空间计算了。

alter table polygons add geom_area  numeric;
update polygons set geom_area = st_area(geom);


create table touching_pairs as
select  a.id as a_id,a.id_array, b.id as b_id from polygons a
inner join polygons b on st_touches(a.geom, b.geom);
create index touching_pairs_id on touching_pairs(a_id);

现在找到接触的多边形。这些 id 触及数组中已有的 id 之一，而该 id 本身不在数组中。将其添加到数组中并将面积相加。

create table pairs as
select distinct on(sort(p.id_array || t.id_array,'asc'))  p.geom_area + a.geom_area as geom_area,
sort(p.id_array || t.id_array,'asc') as id_array
 from polygons p -- only need to change this!
 inner join touching_pairs t
 on t.b_id = any(p.id_array) and t.a_id != all(p.id_array)
 inner join polygons a on a.id = t.a_id
 order by sort(p.id_array || t.id_array,'asc') , t.b_id;

然后您可以使用相同的查询，只是使用对...

 create table trios as
select distinct on(sort(p.id_array || t.id_array,'asc'))  p.geom_area + a.geom_area as geom_area,
sort(p.id_array || t.id_array,'asc') as id_array
 from pairs p -- only need to change this!
 inner join touching_pairs t
 on t.b_id = any(p.id_array) and t.a_id != all(p.id_array)
 inner join polygons a on a.id = t.a_id
 order by sort(p.id_array || t.id_array,'asc') , t.b_id;

然后：

select * from trios where geom_area = 7;

给你;

{2,3,11}
{3,10,11}
{4,10,11}

四重奏...

select distinct on(sort(p.id_array || t.id_array,'asc'))  p.geom_area + a.geom_area as geom_area,
sort(p.id_array || t.id_array,'asc') as id_array
 from trios p -- only need to change this!
 inner join touching_pairs t
 on t.b_id = any(p.id_array) and t.a_id != all(p.id_array)
 inner join polygons a on a.id = t.a_id
 order by sort(p.id_array || t.id_array,'asc') , t.b_id;

我认为这是最快的方法。只需运行它在几秒钟内在 491 个多边形上。

Postgis：查询以获取面积总和等于42的多个相邻多边形的列表

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

postgresql

postgis