1

I would like to perform spatial identity operation (ArcGIS example here) using PostGIS. I have two polygon tables, one of which is the input layer and the other table is the overlay layer. The output should be the original input layer but broken up to pieces where the overlay layer intersects it. This could get complex as multiple overlay polygons can intersect the input polygons. The output table need to have all the original fields of the input table plus a column from the overlay table. This extra column will have the data from the overlay table where they intersect and a fixed value like -1 or NA in non-overlapping areas.

How can this be done in PostGIS?

I am looking for an efficient solution as my input table has about 400k and my overlay table has about 20 million records.

  • ST_Intersects or ST_Intersection would handle this – DPSSpatial Aug 23 at 14:54
1

Install the PostGIS Addons and look at the ST_SplitAgg() example.

Here is a query doing an overlay between two independent tables and joining the attributes of the second layer to the first layer:

WITH geomtableA AS (
  SELECT 0 id, ST_GeomFromText('POLYGON((5 5, 5 7, 7 7, 7 5, 5 5))') geom
  UNION ALL
  SELECT 1 id, ST_GeomFromText('POLYGON((0 0, 0 2, 2 2, 2 0, 0 0), (0.2 0.5, 0.2 1.5, 0.8 1.5, 0.8 0.5, 0.2 0.5))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((2 0.5, 2 1.5, 4 1.5, 4 0.5, 2 0.5))') geom
), geomtableB AS (
  SELECT 3 id, ST_GeomFromText('POLYGON((1 0.2, 1 1, 3 1, 3 0.2, 1 0.2))') geom
  UNION ALL
  SELECT 4 id, ST_GeomFromText('POLYGON((1.5 1, 1.5 1.4, 2.5 1.4, 2.5 1, 1.5 1))') geom
), overlay AS (
  SELECT DISTINCT ON (geom) a.id, unnest(ST_SplitAgg(a.geom, b.geom, 0.00001)) geom
  FROM geomtableA a LEFT JOIN
       geomtableB b
  ON ST_Equals(a.geom, b.geom) OR
        ST_Contains(a.geom, b.geom) OR
        ST_Contains(b.geom, a.geom) OR
        ST_Overlaps(a.geom, b.geom)
  GROUP BY a.id
  ORDER BY geom, max(ST_Area(a.geom)) DESC
)
SELECT a.id aid, b.id bid, a.geom
FROM overlay a LEFT JOIN geomtableB b
ON ST_Within(ST_Centroid(a.geom), b.geom);
  • Thanks @PierreRacine. This is very close to identity operation. So assuming geomtableA as input and geomtableB as overlay, the output above only includes the intersecting split pieces of the two polygons. The output of identity needs to have all the pieces of the input. Now, it needs to add the non-intersecting pieces (with a value like 'null') to the output. The overlay CTE has all the pieces. It should also include all initial non-intersecting polygons in the output with null values (eg POLYGON((5 5, 5 7, 7 7, 7 5, 5 5)). I tried but wasn't able to perform it properly – Kaplan Aug 29 at 5:43
  • also it would be great if you could explain the mechanics of the query. specially why you are not using ST_intersects instead using the four function in the first where clause and ST_Within in the second one. It seems that ST_SplitAgg returns an array of geometries. Are they stacked indentical geometries and hence the need for DISTINCT ON? – Kaplan Aug 29 at 5:50
  • Right! I edited the query with a LEFT JOIN so now it looks more like a identity operation. – Pierre Racine Aug 29 at 20:04
  • Why I don't use ST_Intersects() is a good question. I don't remember why... Was it for processing speed or because there was some problem with overlapping borders? I don't know. Probably you can replace ST_Equals(a.geom, b.geom) OR ST_Contains(a.geom, b.geom) OR ST_Contains(b.geom, a.geom) OR ST_Overlaps(a.geom, b.geom) with ST_Intersects(a.geom, b.geom) without problem. – Pierre Racine Aug 29 at 20:14
  • ST_SplitAgg() splits a polygon with all its intersecting polygons, returning the polygon splitted in many parts. You might ends up with many identical polygons because one splitted polygon part will be identical to another splitted polygon part. That's why you have to DISTINCT ON (geom). The last part of the query just LEFT JOIN each polygon back with the overlaying ones to get their IDs. The beauty of this method is that it works whatever the number of polygons overlapping each other. – Pierre Racine Aug 29 at 20:15
1

These steps should produce an output similiar to the ESRI Identity Tool.
You should create a Spatial Index on each layer, before each command.
I'm unsure if PostGIS creates a spatial index automatically.

  1. ST_Intersection Input and Overlay (this produces stacked polygons that overlap), output as intersect01.
  2. Perform a spatial join on intersect01, adding the 01 Overlay attribute to the 01 Input polygon in intersect01.
  3. Delete the Overlay polygons from intersect01, they are no longer needed.
  4. ST_Difference Input and Intersect01 (these are Input areas that do not overlap Overlay), output as difference02.
  5. Insert difference02 into intersect01, output as Final03.
  • What if two polygons from overlay (o1 and o2) intersects with one polygon from input (i1)? ST_Difference() will remove them independently in two steps (i1 - o1 and i1 - o2). You will end up with overlaps in Final03. – Pierre Racine Aug 23 at 21:04
  • @Pierre, ST_DIfference will erase O from I where they overlap, producing the compliment to I from #1. You should have no extra overlapping I in the Final output. I think you are describing a situation, in the ESRI picture, of 2 circles overlapping the square. This should not be a problem. – klewis Aug 23 at 23:06
  • @Pierre, I've just tested ESRI Identity, it does not produce stacked polygons, the polygons are broken but with all attributes combined. My implementation has stacked polygons from st_Intersection, this is incorrect. – klewis Aug 23 at 23:35

Your Answer

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

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