Cutting Polygons with Polygons in PostGIS – Advanced Techniques

clippolygonpostgispostgresql

I had a really good answer to my question Using CASE to select between two geometry functions? but a problem arises if there are more than one lake in a forest; then the solution returns multiple forests.

Let me rephrase the question in much more general terms.

Problem

Consider this situation
I have some forest polygons, and some lake polygons, in separate tables (layers) in PostgreSQL. I want to use PostGIS/PostgreSQL to create a new table with forests, that have holes where the lakes are.

Some forests have many lakes, some have one, and many have none. I expect the number of output forests to be the same as the input forests (except if a forest is actually cut entirely by a lake).

I need to preserve an ID-number attribute on the forests.

Best Answer

This is essentially the same as the other answers, however if you are dealing with larger tables and cutting lakes out of many forests you may want to try this variation.

It should take advantage of spatial indexes and only union together lakes that need to be used as a cutter.

I've done this as a select (with a CTE to provide sample data)

/* Sample data for forest */
WITH forest AS (
    SELECT *
    FROM (VALUES
        (1,ST_GeomFromText('POLYGON((0 0,10 0,10 10,0 10,0 0))',0))
        ,(2,ST_GeomFromText('POLYGON((10 0,20 0,20 10,10 10,10 0))',0))
        ) Forest(id, geom)
    ),
/* Sample data for lake */
    lake AS (
    SELECT *
    FROM (VALUES
        (1,ST_GeomFromText('POLYGON((2 2,4 2,4 4,2 4,2 2))',0)) /* hole in first */
        ,(2,ST_GeomFromText('POLYGON((8 5,12 5,12 9,8 9,8 5))',0)) /* overlapping */
        ,(3,ST_GeomFromText('POLYGON((12 2,14 2,14 4,12 4,12 2))',0)) /* hole in second */
        ) lake(id, geom)
    )
/* the actual query */
SELECT id, 
    ST_AsText(
        ST_Difference(
            f.geom,
            /* Correlated subquery to fetch only the lakes intersected by the current forest */
            (
                SELECT ST_Union(l.geom) 
                FROM lake l 
                WHERE ST_Intersects(l.geom,f.geom)
            )
        )
    )
FROM forest f
Related Question