[GIS] subtract polygons from other polygons with SQL

spatial-databasesql server

I have a MSSQL table that houses a bunch of polygons. I want to select a few ( with a WHERE clause ) and subtract another bunch from these first ones.

Compare it with having all plot-outlines and building-footprints in a table.

How can i get the difference between these two sets? ( in the example, return all gardens (where a garden = plot – house)

Something along these lines: (this is obviously not working – syntax is wrong too)

SET plots = SELECT geom FROM tableA WHERE type= 'plot'
SET houses = SELECT geom FROM tableA WHERE type= 'house' 

SELECT plots.STDifference(houses) AS gardens FROM tableA

I can't change the data structure, so proposing to move plots and houses into a separate table is not helpful.

Best Answer

If you are working with SQL Server 2012+ then you should be able to do this within a single query. Don't expect it to be very performent, but do make sure that you have all the appropriate indexes (especially spatial) in place.

I'm working on the assumption that you may have multiple hole making polygons that may or may not be entirely within the polygons to be cut.

The following example uses a cross apply to select and union aggregates all the polygons which intersect the target polygons (one at a time). This leaves us with 2 geometries that can be differenced.

The result should be a geometry for each plot polygon with any house polygon cut out of it.

-- Set up a test table
CREATE TABLE tableA (
    id int,
    geom geometry,
    type varchar(5)
    );

INSERT INTO tableA 
VALUES
 (1,Geometry::STGeomFromText('POLYGON(( 0 0, 10 0, 10 10, 0 10, 0 0 ))',0),'plot')
,(2,Geometry::STGeomFromText('POLYGON(( 10 0, 20 0, 20 10, 10 10, 10 0 ))',0),'plot')
,(3,Geometry::STGeomFromText('POLYGON(( 0 10, 10 10, 10 20, 0 20, 0 10 ))',0),'plot')
,(4,Geometry::STGeomFromText('POLYGON(( 2 2, 5 2, 5 5, 2 5, 2 2 ))',0),'house')
,(5,Geometry::STGeomFromText('POLYGON(( 12 2, 15 2, 15 5, 12 5, 12 2 ))',0),'house')
,(6,Geometry::STGeomFromText('POLYGON(( 17 2, 19 2, 19 5, 17 5, 17 2 ))',0),'house')
,(7,Geometry::STGeomFromText('POLYGON(( 2 12, 5 12, 5 15, 2 15, 2 12 ))',0),'house')
,(8,Geometry::STGeomFromText('POLYGON(( 5 9, 7 9, 7 11, 5 11, 5 9 ))',0),'house');

-- Return the plots with houses removed
SELECT p.id
    , p.type
    , p.geom.STDifference(x.exclude) garden -- plot with houses removed
FROM tableA p
    CROSS APPLY ( -- For each plot do the following query
        -- a union aggregate of houses that intersect the current plot
        SELECT Geometry::UnionAggregate(geom) exclude 
        FROM tableA b 
        WHERE type = 'house' 
            AND p.geom.STIntersects(b.geom) = 1
        ) x
WHERE p.type = 'plot'; -- consider plots only for difference

-- Cleanup
DROP TABLE tableA;
Related Question