Spatial Overlay of two Polygon Coverage Tables

Bill Thoen ​asks:

"Using PopstgreSQL/PostGIS, what are the steps I need to use to split one layer of polygons with another layer and also maintain data attrributes among the polygons that are split? For example, I've got corn and soybean field polygons that need to be split at the borders of overlaying CLU (Common Land Unit) polygons, and if a corn field crosses under a CLU boundary, I want to split the field at the CLU boundary and copy the orignal field id, crop code and other attributes into the two (or more) resulting farm field polygons (and probably assign a new unique polygon id to them all too.) "

Paul Ramsey ​answers:

CREATE TABLE new_fields AS
SELECT
   Intersection(f.the_geom, c.the_geom) AS the_geom,
   f.attr1,
   f.attr2,
   c.clu_name
FROM
   fields f,
   clu c
WHERE
   f.the_geom && c.the_geom
AND
   Intersects(f.the_geom, c.the_geom)

Kevin Neufeld answers:

Given two polygonal datasets, fields and clu boundaries,

One could perform an overlay operation by:

1. Extract the linework from the polygons (you may wish to add any interior rings to this query).

CREATE TEMP TABLE all_lines AS
SELECT St_ExteriorRing(the_geom) AS the_geom
FROM fields
UNION ALL
SELECT St_ExteriorRing(the_geom) AS the_geom
FROM clu;

2. Node the linework:

CREATE TEMP TABLE noded_lines AS
SELECT St_Union(the_geom) AS the_geom
FROM all_lines;

3. Re-polygonize all linework:

CREATE TEMP TABLE new_polys (id serial PRIMARY KEY, the_geom geometry);
INSERT INTO new_polys (the_geom)
SELECT geom AS the_geom
FROM St_Dump((
 SELECT St_Polygonize(the_geom) AS the_geom
 FROM noded_lines
));

4. Transfer attributes from original polygons to newly formed polygons
   • Create point-in-polygon features.

CREATE TEMP TABLE new_polys_pip AS
SELECT id, ST_PointOnSurface(the_geom) AS the_geom
FROM new_polys;

• Overlay points with original polygons, transferring the attributes.

CREATE TEMP TABLE pip_with_attributes AS
SELECT a.id, f.attr1, c.attr2
FROM new_polys_pip a
LEFT JOIN fields f ON St_Within(a.the_geom, f.the_geom)
LEFT JOIN clu c ON St_Within(a.the_geom, c.the_geom);

• Join the points to the newly formed polygons, transferring the attributes:

CREATE TABLE new_fields AS
SELECT *
FROM new_polys a LEFT JOIN pip_with_attributes b USING (id);

Overlay Function

Chris Hodgson adds:

Here is a basic SQL function which implements Kevin's multi-step approach in a single query, accessible by one function call:

CREATE OR REPLACE FUNCTION my_overlay(table1 varchar, col1 varchar, table2 varchar, col2 varchar, result_table varchar) returns void AS 
$$ 
BEGIN
EXECUTE 'CREATE TABLE ' || result_table || ' AS '
|| 'WITH all_lines AS ( '
	|| 'SELECT St_ExteriorRing(st_geometryn(the_geom,1)) AS the_geom '
	|| 'FROM ' || table1
	|| ' UNION ALL '
	|| 'SELECT St_ExteriorRing(st_geometryn(the_geom,1)) AS the_geom '
	|| 'FROM ' || table2
|| '), noded_lines AS ( '
	|| 'SELECT St_Union(the_geom) AS the_geom '
	|| 'FROM all_lines '
|| '), new_polys AS (  '
	|| 'SELECT geom AS the_geom, ST_PointOnSurface(geom) AS pip '
	|| 'FROM St_Dump(( '
		|| 'SELECT St_Polygonize(the_geom) AS the_geom '
		|| 'FROM noded_lines)) '
|| ') '
|| 'SELECT a.' || col1 || ' as ' || table1 || '_' || col1 
|| ', b.' || col2 || ' as ' || table2 || '_' || col2 
|| ', p.the_geom as the_geom '
|| 'FROM new_polys p '
|| 'LEFT JOIN ' || table1 || ' a ON St_Within(p.pip, a.the_geom) '
|| 'LEFT JOIN ' || table2 || ' b ON St_Within(p.pip, b.the_geom)';
END
$$ LANGUAGE plpgsql;

The function takes the name of two tables, and the names of a column from each of them (the id/primary key is the expected use) and the new table name to be created. I believe this approach will be memory limited because all of the linework, and the resulting polygons, are in memory at once and stored in a single (potentially very large) multi-geometry.

Note that the result set will include resultants for areas which are in holes but which do not overlap any parent polygon. If not desired these can be filtered out by keeping only records for which at least one id is non-NULL.

Worked Example

This is the SQL from the above function, with example data:

WITH poly_a(id, geom) AS (VALUES
    ( 'a1', 'POLYGON ((10 40, 30 40, 30 10, 10 10, 10 40))'::geometry ),
    ( 'a2', 'POLYGON ((70 10, 30 10, 30 90, 70 90, 70 10), (40 40, 60 40, 60 20, 40 20, 40 40), (40 80, 60 80, 60 60, 40 60, 40 80))'::geometry ),
    ( 'a3', 'POLYGON ((40 40, 60 40, 60 20, 40 20, 40 40))'::geometry )
)
,poly_b(id, geom) AS (VALUES
    ( 'b1', 'POLYGON ((90 70, 90 50, 50 50, 50 70, 90 70))'::geometry ),
    ( 'b2', 'POLYGON ((90 30, 50 30, 50 50, 90 50, 90 30))'::geometry ),
    ( 'b2', 'POLYGON ((90 10, 70 10, 70 30, 90 30, 90 10))'::geometry )
)
,lines AS ( 
  SELECT ST_Boundary(geom) AS geom FROM poly_a
  UNION ALL
  SELECT ST_Boundary(geom) AS geom FROM poly_b
)
,noded_lines AS ( SELECT St_Union(geom) AS geom FROM lines ) 
,resultants AS (  
  SELECT geom, ST_PointOnSurface(geom) AS pip 
    FROM St_Dump(
           ( SELECT ST_Polygonize(geom) AS geom FROM noded_lines ))   
)
SELECT a.id AS ida, b.id AS idb, r.geom
  FROM resultants r
  LEFT JOIN poly_a a ON St_Within(r.pip, a.geom) 
  LEFT JOIN poly_b b ON St_Within(r.pip, b.geom)
  WHERE a.id IS NOT NULL OR b.id IS NOT NULL;