Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a PostGIS table of polygons where some intersect with one another. This is what I'm trying to do:

  • For a given polygon selected by id, give me all of the polygons that intersect. Basically, select the_geom from the_table where ST_Intersects(the_geom, (select the_geom from the_table where source_id = '123'))
  • From these polygons, I need to create a new polygons such that intersection becomes a new polygon. So if polygon A intersects with polygon B, I will get 3 new polygons: A minus AB, AB, and B minus AB.

Any ideas?

share|improve this question
1  
Hmmm, I think see what you're getting at but a simple graphic might do wonders to help me (and others) see exactly what you want. –  Jason Jul 23 '10 at 2:19
    
Added some in my answer. –  Adam Matan Aug 2 '10 at 7:37
add comment

4 Answers

up vote 10 down vote accepted

Since you said you get a group of intersecting polygons for each polygon you're interested in, you may want to create what is referred to as a "polygon overlay".

This isn't exactly what Adam's solution is doing. To see the difference, take a look at this picture of an ABC intersection:

ABC intersection

I believe Adam's solution will create an "AB" polygon that covers both the area of "AB!C" and "ABC", as well as an "AC" polygon that covers "AC!B" and "ABC", and a "BC" polygon that is "BC!A" and "ABC". So the "AB", "AC", and "BC" output polygons would all overlap the "ABC" area.

A polygon overlay produces non-overlapping polygons, so AB!C would be one polygon and ABC would be one polygon.

Creating a polygon overlay in PostGIS is actually pretty straightforward.

There are basically three steps.

Step 1 is extract the linework [Note that I'm using the exterior ring of the polygon, it does get a little more complicated if you want to correctly handle holes]:

SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines

Step 2 is to "node" the linework (produce a node at every intersection). Some libraries like JTS have "Noder" classes you can use to do this, but in PostGIS the ST_Union function does it for you:

SELECT ST_Union(the_geom) AS the_geom FROM (...your lines...) AS noded_lines

Step 3 is to create all the possible non-overlapping polygons that can come from all those lines, done by the ST_Polygonize function:

SELECT ST_Polygonize(the_geom) AS the_geom FROM (...your noded lines...)

You could save the output of each of those steps into a temp table, or you can combine them all into a single statement:

CREATE TABLE my_poly_overlay AS
SELECT geom FROM ST_Dump((
    SELECT ST_Polygonize(the_geom) AS the_geom FROM (
        SELECT ST_Union(the_geom) AS the_geom FROM (
            SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines
        ) AS noded_lines
    )
)

I'm using ST_Dump because the output of ST_Polygonize is a geometry collection, and it is (usually) more convenient to have a table where each row is one of the polygons that makes up the polygon overlay.

share|improve this answer
add comment

If I understand correctly, You want to take (A is the left geometry, B is the right):

Image of A∪B

And extract:

A∖AB

Image of A∖AB

AB

Image of AB

and B∖AB

Image of B∖AB

That is - three different geometries for every intersecting pair.

First, let's create a view of all intersecting geometries. Assuming your table name is polygons_table, we will use:

CREATE OR REPLACE VIEW p_intersections AS    -- Create a view with the 
SELECT t1.the_geom as t1_geom,               -- intersecting geoms. Each pair
       t2.the_geom as t2_geom                -- appears once (t2.id<t2.id)
    FROM polygons_table t1, polygons_table t2  
         WHERE t1.id<t2.id AND t1.the_geom && t2.the_geom 
                           AND intersects t1.the_geom, t2.the_geom;

Now we have a view (practically, a read-only table) that stores pairs of intersecting geoms, where each pair appears only once due to the t1.id<t2.id condition.

Now let's gather your intersections - A∖AB,AB and B∖AB, using SQL's UNION on all three queries:

--AB:
SELECT ST_intersection(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION 

--A∖AB:
SELECT ST_Difference(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION

--B∖AB:
SELECT ST_Difference(t2.the_geom, t1.the_geom) 
    AS geom 
    FROM p_intersections;

Notes:

  1. The && operator is used as a filter before the intersects operator, to improve performance.
  2. I've chosen to create a VIEW instead of one gigantic query; This is for convenience only.
  3. If you meant AB is the union, not the intersection, of A and B - Use ST_Union instead of st_intersection at the UNION query in the appropriate places.
  4. The sign is a unicode sign for Set difference; remove it from the code if it confuses your database.
  5. Pictures courtesy of Wikimedia Commons' nice Set theory category.
share|improve this answer
    
My bug ticket on meta: meta.gis.stackexchange.com/questions/79/… –  Adam Matan Aug 2 '10 at 7:50
    
Nice explanation! Results are the same as in scw solution, but his way should be faster (does not compute /or store/ additional intersections of A and B) –  stachu Aug 3 '10 at 6:51
    
Thanks! I think I don't store any extra information, as I only create SQL VIEWs, not tables. –  Adam Matan Aug 3 '10 at 6:56
    
Yes, that is true, but you compute additional Intersection of A and B, which is not necessery –  stachu Aug 3 '10 at 12:21
1  
Images in this answer don't work anymore. –  Fezter Feb 2 '13 at 6:39
show 2 more comments

What you're describing is the way that the Union operator works in ArcGIS, but its a little different than either Union or Intersection in the GEOS world. The Shapely manual has examples of how sets work in GEOS. However, the PostGIS wiki does have a good example using linework which should do the trick for you.

Alternatively, you could compute three things:

  1. ST_Intersection(geom_a, geom_b)
  2. ST_Difference(geom_a, geom_b)
  3. ST_Difference(geom_b, geom_a)

Those should be the three polygons you mentioned in your second bullet point.

share|improve this answer
    
The PostGIS wiki example is good –  fmark Jul 23 '10 at 2:41
    
Wouldn't ST_Intersects return boolean if they intersect or not? I think ST_Intersection would work. –  Jason Jul 23 '10 at 2:52
    
Yeah, typo on my part -- fixed in the original now, thanks Jason! –  scw Jul 23 '10 at 2:56
add comment

Something like:

INSERT INTO new_table VALUES((select id, the_geom from old_table where st_intersects(the_geom,(select the_geom from old_table where id='123')) = true

EDIT: you need the actual intersection of the polygon.

INSERT INTO new_table values((select id, ST_Intersection(the_geom,(select the_geom from old where id = 123))

see if that works out.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.