3

I am attempting to replicate the ArcGIS erase tool in PostGIS.

Here is the problem I am working with

Two Layers

  1. land use land cover (150,000) multipolygon table
  2. multipolygon soils table (30,000)

    • both have spatial indexes

here is my query

CREATE TABLE agriculture.important_ag_soils_undeveloped AS
--cte selecting farmland data
with farm_to_clip AS (
  SELECT * 
  FROM base_layers.ssurgo_soils_highlands_updated_stateplane 
  WHERE farmlndcl IN ('Farmland of statewide importance','Farmland of unique importance','All areas are prime farmland')
),
--cte selecting land use land cover data to erase from farm_to_clip data above
lulc_to_erase AS (
  SELECT shape FROM base_layers.lulc_huc_hl_clip 
  WHERE lu12 IN (1110,1120,1130,1140,1150,1200,1300,1500,1600,7300,1211,1400,1410,1411,1419,1420,1440,1810,1462,1463,1499,1710,1800,1804))
--now I am running this sort of messy query where I perform a st_difference for each farm polygon record and run a subquery as the second input unioning the landuse land cover where it intersects the farm polygon 
-- which yeilds the difference, and I wrap the st_difference function in the coalesce function so if it is st_difference(f.shape,NULL) then f.shape is returned for that row. 
SELECT 
  f.*,
  coalesce(
    st_difference(f.shape,
      (SELECT st_union(l.shape) 
       FROM lulc_to_erase AS l 
       WHERE st_intersects(f.shape,l.shape)
      )
    ),
    f.shape) AS geom
FROM farm_to_clip f;

this query takes about 5 minutes to run and I know there must a more efficient way of writing this query.

from @tilt suggestion I removed a the subselect and coalesce function and just used a left join with the st_difference function here is my updated query:

create table agriculture.important_ag_soils_undeveloped2 as

WITH farm_to_clip AS (
    SELECT * 
    FROM base_layers.ssurgo_soils_highlands_updated_stateplane
    WHERE farmlndcl IN ('Farmland of statewide importance','Farmland of unique importance','All areas are prime farmland')
    ),

lulc_to_erase AS (
    SELECT shape from base_layers.lulc_huc_hl_clip 
    WHERE lu12 in (1110,1120,1130,1140,1150,1200,1300,1500,1600,7300,1211,1400,1410,1411,1419,1420,1440,1810,1462,1463,1499,1710,1800,1804))

SELECT f.*,
    ST_Difference(f.shape,l.shape) geom 
FROM farm_to_clip f 
LEFT JOIN lulc_to_erase l ON ST_Intersects(f.shape,l.shape);

Problem: the new method of using the left join returns a significant more amount of polygons in the new layer than the first query I used.

My question becomes which query gives me the farm layer erased by the lulc like the way ArcGIS performs its erase?

UPDATE

I used the arcgis erase tool on the farm_to_clip and lulc_to_erase and the layer that was created is the same layer that I created in my first query. My second query did not seem to work correctly at all, but I cannot figure out why..?

I am going to try unioning my lulc_to_erase geometry before I send it through the select st_difference function

farm = light transparent green
lulc = hatched orange 

enter image description here

QUERY 1 RESULTS

enter image description here

QUERY 2 RESULTS

enter image description here

  • 1
    Why did you make this subselect inside your ST_Difference? It seems like this could be done with a normal LEFT JOIN and GROUP BY. – tilt Aug 1 '17 at 20:14
  • oh yeah let me try that...and I would still need the coelesce function in case the st_difference's second input (land cover) is null – ziggy Aug 1 '17 at 20:17
  • @tilt I ran a left join but I returned a lot more polygon area than my first query.. now I am unsure which one is more accurate – ziggy Aug 2 '17 at 13:53
  • You should get a feeling of what happened when you visually inspect the output in QGIS or ArcGis. Maybe you can share some pictures of the situation before and after? – tilt Aug 2 '17 at 15:20
  • @tilt updated, the first query seems to work but I cannot wrap my head around why the second query with the left join somehow returns more than the first query. I get that it returns all the farm rows with the left join, but if it a farm poly does not intersect a lulc poly then it will send a null geom into the st_difference function...? – ziggy Aug 2 '17 at 15:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.