3

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

8
  • 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
    – ziggy
    Oct 29, 2016 at 16:47
  • 1
    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
    – Matte
    Oct 30, 2016 at 6:31
  • 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 ...
    – Kurt
    Oct 30, 2016 at 6:37
  • I'd run SELECT DISTINCT ST_GeometryType(st_collect) FROM layer_3; to find out what you're actually dealing with. Oct 30, 2016 at 18:33
  • 1
    @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). Oct 31, 2016 at 6:29

2 Answers 2

5

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.

4
  • 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?
    – Kurt
    Oct 29, 2016 at 20:07
  • 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
    – raphael
    Oct 29, 2016 at 20:48
  • ok, the option is enabled by default. I will run ANALYZE layername on each and proceed
    – Kurt
    Oct 30, 2016 at 4:58
  • according to your instructions I created a smaller test-suite, please see the edited question above. any hints??
    – Kurt
    Oct 30, 2016 at 6:06
0
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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