多边形之间的PostGIS递归相交 [英] PostGIS recursive intersection between polygons

问题描述

我正在尝试在空间表中的所有多边形之间执行递归相交，以获取生成的(多)多边形和每个相交信息.

I'm trying to perform a recursive intersection between all the Polygons in a spatial Table, and obtain the resulting (multi)pololygons and the information about every intersection for each of them.

一张图片(不是真正的比例尺)来解释它:

An image (not really in scale) to explain it:

让我们说一个表中有A, B, C个正方形.我想在输出中包含A, B, C, A+B, A+C, B+C, A+B+C多边形，我需要知道A+B是A和B的交集，依此类推.

Let's say there are A, B, C squares in a table. I'd like to have A, B, C, A+B, A+C, B+C, A+B+C polygons in output, and I need to know that A+B is the intersection of A and B and so on.

到目前为止，我有一个执行相交的查询，但是它没有切除"原始多边形的相交部分.例如:

So far I have a query which performs the intersections, but it does not "cut off" the intersected part of the original polygons. For example:

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多边形结果的图像:

An image of the result I get now for the A and A+C polygons:

这是一个测试脚本，使用图像中的正方形作为数据.查看area列，很明显缺少一些递归的ST_Difference，我只是想不通.任何想法都欢迎.

Here is a test script, using the squares in the images as data. Looking at the area column, it is clear some recursive ST_Difference is missing, I just can't figure out how. Any idea is welcomed.

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

在此特定示例中，窗口函数并不是严格必需的，但是此代码是我的真实案例的简化版本，它在侧面还做了一些其他事情.

The window function isn't strictly necessary in this particular example of course, but this code is a simplified version of my real case, which does a few more things on the side.

我正在使用PostgreSQL 9.2和PostGIS 2.0

I'm using PostgreSQL 9.2 and PostGIS 2.0

推荐答案

ST_DIFFRENCE不必是递归的，您已经拥有了所有多边形，因此对于每个几何，您都必须减去包含该ret的其他几何的并集，不等于它.这行得通，所以您应该像这样:

ST_DIFFRENCE doesn't have to be recursive, you already have all the polygons so from every geom you have to substract the union of other geoms which contain that ret but ain't equal to it. This works so you should do it kinda like that:

    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

这篇关于多边形之间的PostGIS递归相交的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！