11

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.

3
  • 1
    FYI, in ArcGIS this operation is called "Erase", your blue layer is "erasing" the green one. In QGIS, I think it's called "Difference". This may help: gis.stackexchange.com/questions/90174/…
    – Dan C
    Oct 11, 2016 at 19:30
  • 1
    Yes, the equivalent function in Esri is 'erase'. And the solution in QGIS likely involves ST_Difference(). The link I refer to also use ST_Difference(). But I need to extend it a bit more to make it work through an entire table. And the 'tricker' based solution is not applicable, as I'm not manually editing, but simply trying to Esri-erase one layer with another.
    – Martin
    Oct 11, 2016 at 20:00
  • 1
    Did you try a ST_UNION on the geometry that should be subtracted (in your case, the lakes)?
    – Greg Z
    Oct 11, 2016 at 20:03

5 Answers 5

10

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
4

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

3
  • Shouldnt it be FROM forest, union_lakes as last line?
    – BERA
    Nov 13, 2018 at 12:49
  • Not if you only want the forests returned. The union_lakes table is used in the ST_Difference, but otherwise ignored for the returned data. Nov 13, 2018 at 20:50
  • 1
    Ok Im getting an error like "missing from clause", but I could be doing something wrong
    – BERA
    Nov 14, 2018 at 10:04
2

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

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])
);
1

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... :-)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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