[GIS] How to intersect two polygon-layers in PostgreSQL/PostGIS

intersectionpostgispostgresql

EDIT:

according to the instructions given by @raphael I created a new, smaller test-set. the old, obsolete original question is now at the bottom of this post

1) Imported two shapefiles into postgresql/postgis:

layer1 with 39 polygons

layer2 with 79 polygons

I used the pgadminIII plugin for the import of the shapefiles. Option is enabled for automatically generating a spatial index.

2) ran

analyze layer1;
analyze layer2;

3) ran

select *,ST_IsValidReason(geom)  from layer1
where ST_IsValid(geom) = false;
select *,ST_IsValidReason(geom)  from layer2
where ST_IsValid(geom) = false;

both layers are valid!

now I have 2 tables with spatialtype MultiPolygon in the database

see screenshot2:
layer1 = blue
layer2 = yellow
enter image description here
4) now ran

create table layer3 
as
SELECT layer1.id, layer1.bev, 
       ST_COLLECT(ST_INTERSECTION(layer1.geom, layer2.geom))
FROM layer1, layer2
WHERE ST_Intersects(layer1.geom, layer2.geom)
GROUP BY layer1.id,layer1.bev;

5) ran

select sum(bev) from layer1;

returns 56462

and 

select sum(bev) from layer3;

also returns 56462

select count(*) from layer3;

returns 39

so all seems ok

but when I import layer3 into qgis I only have 27 polygons !!
(see screenshot 3)
enter image description here

6) ran

select *,ST_IsValidReason(st_collect)  from layer3
where ST_IsValid(st_collect) = false;

delivers no invalid rows !!!

and

select count(st_collect) from layer3;

delivers 39

so what I am doing wrong?

edit2:

the following seems to work:

create table layer4
as
SELECT layer1.id, layer1.bev, 
      st_union(ST_INTERSECTION(layer1.geom, layer2.geom))
FROM layer1, layer2
WHERE ST_Intersects(layer1.geom, layer2.geom)
GROUP BY layer1.id,layer1.bev;

but I must confirm, that I dont really understand why …, a simple explanation for a mere mortal ??

old, now obsolete original question:
I have 2 layers:

layer #1 named burgenlandzsp has 382 polygons (= red in screenshot)

layer #2 named siedlungsraum_singleparts_burgenland has 897 polygons (= green in screenshot)

Problem: I need to create a new layer which contains only the parts of layer#1 which intersect with layer #2.

I tried to do this in PostGIS with the following command:

create table bglsiedlungsraum as
select
   burgenlandzsp.id,burgenlandzsp.name,burgenlandzsp.zsp,
   burgenlandzsp.bev,
   ST_intersection(burgenlandzsp.geom,
                   siedlungsraum_singleparts_burgenland.geom) as geom 
from burgenlandzsp,siedlungsraum_singleparts_burgenland
where ST_intersects(burgenlandzsp.geom, 
                    siedlungsraum_singleparts_burgenland.geom);

I expected to get a new table(layer) with a multigeometry of 382 multipolygons, but it did not work

after about one hour this query "crashed" my hp-laptop (6mb ram, win10 64 bit, 4 cores @1,6 GHz)

enter image description here

Best Answer

I expected to get a new table(layer) with a multigeometry of 382 multipolygons, but it did not work

Your expectation is incorrect. PostGIS is only returning ANY combination of layer_1 and layer_2 that intersects. If every polygon in layer_2 straddled two polygons in layer_1 the result of 

SELECT COUNT(*)
FROM layer_1, layer_2
WHERE ST_Intersects(layer_1.geom, layer_2.geom)

Would be 897 * 2.

If you want 382 multipolygons you need to use a combination of ST_Collect with GROUP BY like so

SELECT layer_1.id, ST_COLLECT(ST_INTERSECTION(layer_1.geom, layer_2.geom))
FROM layer_1, layer_2
WHERE ST_Intersects(layer_1.geom, layer_2.geom)
GROUP BY layer_1.id

after about one hour this query "crashed" my hp-laptop

Are you using a spatial index on either layer? If not then PostGIS is doing 382 * 897 comparisons to check if each layer_1 polygon intersects with each layer_2 polygon, which may take a while.

Related Question