More on cutting polygons with polygons, in PostGIS

Question

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

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.

Answer 1

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

Answer 2

As mentioned by Greg Z in a comment, you first need to union your lakes so that they become one geometry (of type MULTIPOLYGON). I would do with with a common table expression (CTE):

WITH union_lakes AS (
    SELECT ST_Union(lake.geom) as geom
    FROM lake
    )
SELECT forest.*, ST_Difference(forest.geom, union_lakes.geom) 
FROM forest

I can't test this query right now, but I believe that should work. If this leads to poor performance, it might be worthwhile taking the CTE out and storing the result as its own table (with its own spatial index).

Answer 3

I just had to clip the ZIP codes that fall within the Denver County boundary and modified the SQL on the ST_Intersection documentation - the results worked perfectly!

select z.geoid10
, ST_Buffer(
ST_Intersection(c.geom, z.geom)
, 0.0)
as geom

from 
ccd."CO_Counties" as c
inner join ccd."CO_ZipCodes" as z
on ST_Intersects(c.geom, z.geom)

where c.name = 'Denver'

Answer 4

Just tested this part, ST_Difference and ST_Union should do the work:

UPDATE forest SET geom_forest = ST_Difference (
   geom_poly, 
   (SELECT ST_UNION(geom_lakes) FROM lakes)
);

with

CREATE TABLE forest
(
  id serial PRIMARY KEY,
  geom_forest geometry(Polygon,[YOUR SRID])
);

CREATE TABLE lakes
(
  id serial PRIMARY KEY,
  geom_lakes geometry(Polygon,[YOUR SRID])
);

Answer 5

I have used your replies, and pieced together some code that does exactly what I was looking fore. It looks like this:

DROP TABLE IF EXISTS samp.lake_cuts_into_forest;

-- Get all the forest with lakes, with the lakes cut out
CREATE TABLE samp.lake_cuts_into_forest AS
  WITH
    forest AS (SELECT * FROM samp.gdk_skov_deep),
    lake AS (SELECT * FROM samp.gdk_soe_deep)
  SELECT DISTINCT forest.gmlid AS forest_gmlid,
    ST_Difference(
      forest.wkb_geometry,
      ( SELECT ST_Union(lake.wkb_geometry) -- Correlated subquery to fetch only the lakes intersected by the current forest
        FROM lake
        WHERE ST_Intersects(lake.wkb_geometry,forest.wkb_geometry)
      )
    ) AS geom
FROM forest, lake
WHERE ST_Intersects(lake.wkb_geometry,forest.wkb_geometry); -- otherwise give blank geom columns for the no-hits

-- Now add all the forest that have no lakes
WITH
    forest AS (SELECT * FROM samp.gdk_skov_deep),
    lake AS (SELECT * FROM samp.gdk_soe_deep),
    all_lakes AS (SELECT ST_Union(lake.wkb_geometry) AS geom FROM lake
)
INSERT INTO samp.lake_cuts_into_forest
  SELECT forest.gmlid AS forest_gmlid,
    forest.wkb_geometry as geom
  FROM forest, all_lakes
  WHERE NOT ST_Intersects(all_lakes.geom,forest.wkb_geometry)

Thanks to all... :-)