How to intersect two polygon-layers in PostgreSQL/PostGIS?

Question

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

• ziggy: im assuming your bringing this layer into QGIS? maybe try doing a select into statement instead of a create table as..that has worked for me before but I cannot tell you why
• Matte: Could be a problem of geometry types. Could be that postgis delivers Multipolygons and Qgis assumes Polygon. The function ST_Dump will expand the Multipolygons to Polygons. Usage is: ST_Dump(actual geometry) as new_geometry
• Kurt: ok, trying to export via the pgadminIII plugin to a shapefile I get the error-message : incompatible mixed geometry types in the table. will look into the documentation for st_dump ...
• Michal Zimmermann: I'd run SELECT DISTINCT ST_GeometryType(st_collect) FROM layer_3; to find out what you're actually dealing with.
• Michal Zimmermann: @Kurt go through the docs, the main difference between those two is that ST_Collect won't dissolve the feature boundaries while ST_Union will (or at least will give it a try and fail due to invalid input).

Answer 1

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.

• Kurt: so the shape-file-import plugin form pgadmin III does not automatically create an index? using the link you thankfully provided, I can create a index : CREATE INDEX [indexname] ON [tablename] USING GIST ( [geometryfield] ); but how do I tell postgis to use this index? how must I modify your example?
• raphael: I think there's an option to create an index, but it's not enabled by default. You need to add the index to both tables, and then run ANALYZE layer_1; on each. Then PostGIS should just use the index automagically
• Kurt: ok, the option is enabled by default. I will run ANALYZE layername on each and proceed
• Kurt: according to your instructions I created a smaller test-suite, please see the edited question above. any hints??

Answer 2

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;

quote from the docs: ST_Collect and ST_Union are often interchangeable. ST_Union is in general orders of magnitude slower than ST_Collect because it tries to dissolve boundaries and reorder geometries to ensure that a constructed Multi* doesn't have intersecting regions.