PostGIS在多边形之间递归交集
我试图在空间表中的所有多边形之间执行递归交集,并获取结果(多)多边形以及每个交点的信息。
一个图像(不是真正的比例尺)来解释它:
假设一张桌子上有A, B, C
格子。 我想在输出中有A, B, C, A+B, A+C, B+C, A+B+C
多边形,并且我需要知道A+B
是A
和B
的交点,等等。
到目前为止,我有一个执行交叉点的查询,但它不会“截断”原始多边形的相交部分。 例如:
Polygon A should be A - (A+B) - (A+C) - (A+B+C)
Polygon A+C should be A+C - (A+B+C)
我现在得到的A
和A+C
多边形的结果图像:
这是一个测试脚本,使用图像中的正方形作为数据。 看看area
列,很明显,一些递归ST_Difference缺失,我只是不知道如何。 任何想法都欢迎。
-- Create a test table
CREATE TABLE test (
name text PRIMARY KEY,
geom geometry(POLYGON)
);
-- Insert test data
INSERT INTO test (name, geom) VALUES
('A', ST_GeomFromText('POLYGON((1 2, 1 6, 5 6, 5 2, 1 2))')),
('B', ST_GeomFromText('POLYGON((0 0, 0 4, 4 4, 4 0, 0 0))')),
('C', ST_GeomFromText('POLYGON((2 0, 2 4, 6 4, 6 0, 2 0))'));
-- Query
WITH RECURSIVE
source (rownum, geom, ret) AS (
SELECT row_number() OVER (ORDER BY name ASC), ST_Multi(geom), ARRAY[name] FROM test
),
r (rownum, geom, ret, incroci) AS (
SELECT rownum, geom, ret, 0 FROM source
UNION ALL
SELECT s.rownum, ST_CollectionExtract(ST_Intersection(s.geom, r.geom), 3), (r.ret || s.ret), (r.incroci + 1)
FROM source AS s INNER JOIN r ON s.rownum > r.rownum AND ST_Intersects(s.geom, r.geom) AND ST_Area(ST_Intersection(s.geom, r.geom)) > 0.5
),
result (geom, ret) AS (
SELECT ST_Union(geom) AS geom, ret FROM r GROUP BY ret
)
SELECT geom, ST_Area(geom) AS area, ret FROM result ORDER BY ret
在这个特定的例子中,窗口函数并不是绝对必要的,但是这个代码是我真实情况的简化版本,它可以做更多的事情。
我使用的是PostgreSQL 9.2和PostGIS 2.0
ST_DIFFRENCE不一定是递归的,你已经拥有所有的多边形,所以你必须减去包含该ret的其他几何体的并集,但不等于它。 这有效,所以你应该这样做:
WITH RECURSIVE
source (rownum, geom, ret) AS (
SELECT row_number() OVER (ORDER BY name ASC), ST_Multi(geom), ARRAY[name] FROM test
),
r (rownum, geom, ret, incroci) AS (
SELECT rownum, geom, ret, 0 FROM source
UNION ALL
SELECT s.rownum, ST_CollectionExtract(ST_Intersection(s.geom, r.geom), 3), (r.ret || s.ret), (r.incroci + 1)
FROM source AS s INNER JOIN r ON s.rownum > r.rownum AND ST_Intersects(s.geom, r.geom) AND ST_Area(ST_Intersection(s.geom, r.geom)) > 0.5
),
result (geom, ret) AS (
SELECT ST_Difference(ST_Union(r.geom),q.geom) AS geom, r.ret FROM r JOIN (SELECT r.ret,ST_UNION(COALESCE(r2.geom,ST_GeomFromText('POLYGON EMPTY'))) as geom FROM r LEFT JOIN r AS r2 ON r.ret<@r2.ret AND r.ret!=r2.ret GROUP BY r.ret) AS q on r.ret=q.ret GROUP BY r.ret,q.geom
)
SELECT geom, ST_Area(geom) AS area, ret FROM result ORDER BY ret
链接地址: http://www.djcxy.com/p/70271.html