4

I am new to postgis and I need to find intersection of three polygons. In my table named "mytable" contains large no of polygon coordinates and grouped into three categories by adding a field called "filename".I labelled the category names as part1, part2, part3 in filename . Now I need to find all polygon which spacially intersect in three categories. I write my code as:

SELECT ST_Intersection(a.geom, b.geom) as intersect_ab
FROM mytable a
INNER JOIN mytable b
ON ST_Intersects(a.geom,b.geom)
INNER JOIN mytable c
ON ST_Intersects(a.geom,c.geom)
WHERE ST_Overlaps(a.geom, b.geom)
AND ST_Overlaps(b.geom, c.geom)
AND ST_Overlaps(a.geom, c.geom)
AND ST_isvalid(a.geom)='t'
AND ST_isvalid(b.geom)='t'
AND ST_isvalid(c.geom)='t'
AND a.filename = 'part1'
AND b.filename = 'part2'
AND c.filename = 'part3');

But it returns wrong answer

4

It sounds like you should probably split your data into three tables first. But if there's a good reason not to do that, then you can use a WITH statement to collect them into table-like entities. If you don't mind me saying, your existing query is quite confusing, so I'm not 100% clear what you are doing. Perhaps if you explain why you want to do this, it would help people answer more thoroughly?

Still, given your existing structure, I think you want a query like this:

WITH a AS (
    SELECT * FROM mytable
    WHERE filename = 'part1'
    AND ST_IsValid(geom::geometry)
    LIMIT 1000
), b AS (
    SELECT * FROM mytable
    WHERE filename = 'part2'
    AND ST_IsValid(geom::geometry)
    LIMIT 1000
), c AS (
    SELECT * FROM mytable
    WHERE filename = 'part3'
    AND ST_IsValid(geom::geometry)
    LIMIT 1000
)
SELECT a.*, b.*, c.*
FROM a, b, c
WHERE ST_Intersects(a.geom, c.geom)
AND ST_Intersects(a.geom, b.geom)
AND ST_Intersects(b.geom, c.geom);

This gets much more complicated with each additional category, and there are certainly better ways to aggregate features like this. This is also a surefire way to issue an incredibly expensive query, so I've put in LIMIT 1000 statements that you can remove when you are ready for the full query.

  • No, It also returns wrong answer.... – A S Feb 20 '17 at 11:18
  • You querstion isn't very clear. How is it the wrong answer? By "Now I need to find all polygon which spacially intersect in three categories" do you mean that the you are only interested in intersections where all 3 categories overlap? – John Powell Feb 20 '17 at 12:10
  • yes. I need only the region which intersect in 3 categories. Assume we have three clusters 1,2 and 3. I need to find regions which are common in all three clusters. – A S Feb 20 '17 at 12:16
  • For example: POLYGON((0 0,0 1,1 1,1 0,0 0)) filename as part1..... POLYGON((0 0,0 3,3 3,3 0,0 0)) filename as part1...... POLYGON((0 0,0 2,2 2,2 0,0 0)) filename as part2...... POLYGON((0 0,0 5,5 5,5 0,0 0)) filename as part2...... POLYGON((0 0,0 3,3 3,3 0,0 0)) filename as part3...... POLYGON((0 0,0 4,4 4,4 0,0 0)) filename as part3....... Now I need to find regions which are common in all categories.***** that is I expect the result as: POLYGON((0 0,0 1,1 1,1 0,0 0)).... POLYGON((0 0,0 2,2 2,2 0,0 0))..... POLYGON((0 0,0 3,3 3,3 0,0 0)) – A S Feb 20 '17 at 12:34
  • Ok I think I understand, I'll edit this shortly. – Richard Law Feb 20 '17 at 18:13
1

Building on what @Richard Law posted earlier, along with the ideas here you should be able to get the intersecting polygons with the following. As with his answer, I'm not sure that this is the most efficient way to accomplish the task.

Since you say you are new to PostGIS, I'll throw out a quick reminder that you want to put a spatial index on your geometry column to improve performance.

WITH a AS (
    SELECT * FROM mytable
    WHERE filename = 'part1'
    AND ST_IsValid(geom::geometry)
    LIMIT 1000
), b AS (
    SELECT * FROM mytable
    WHERE filename = 'part2'
    AND ST_IsValid(geom::geometry)
    LIMIT 1000
), c AS (
    SELECT * FROM mytable
    WHERE filename = 'part3'
    AND ST_IsValid(geom::geometry)
    LIMIT 1000
)

SELECT
    CASE
    WHEN ST_CoveredBy(ab.geom, c.geom)
        THEN ab.geom
    WHEN ST_CoveredBy(c.geom, ab.geom)
        THEN c.geom
    ELSE
        ST_Multi(ST_Intersection(ab.geom, c.geom))
    END as geom
    FROM c,   
        (SELECT 
            CASE
            WHEN ST_CoveredBy(a.geom, b.geom)
                THEN a.geom
            WHEN ST_CoveredBy(b.geom, a.geom)
                THEN b.geom
            ELSE
                ST_Multi(ST_Intersection(a.geom, b.geom))
            END AS geom
        FROM a, b
        ) AS ab;

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.