I am attempting to replicate the ArcGIS erase tool in PostGIS.
Here is the problem I am working with
- land use land cover (150,000) multipolygon table
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?
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
QUERY 1 RESULTS
QUERY 2 RESULTS