Extracting data from four intersecting polygon layers in PostGIS?

Question

I've been using the following code in PostgreSQL/PostGIS for the last year, as needed, which is about once a quarter. It works OK, it's just that it takes days to run. With over 25,000 parcels making up the county, maybe that's fine, but I'd like some expert input to see if there's a better way to do it.

An additional consideration, wetlandyn is a multi-polygon, the rest are single, though soils and luc could be converted to a multi-polygon.

Is querying against multi-polygons any more or less efficient than singles?

SELECT
    parcels.pin AS "PIN",
    soils.musym AS "Soils",
    parcels.camacls as "CAMACLS",
    luc.landcode as "Land Code",
    wetlandyn.inwetland as "In Wetland (y/n)",
    ROUND(((SUM(ST_AREA(ST_INTERSECTION((ST_INTERSECTION(parcels.geom,luc.geom)),
    (ST_INTERSECTION(wetlandyn.geom,soils.geom)))))/ST_AREA(parcels.geom))*parcels.acres)::numeric,4) AS "Legal Soils"

FROM
    soils,
    luc,
    parcels,
    wetlandyn

WHERE
    ST_INTERSECTS(parcels.geom, soils.geom)
    AND
    ST_INTERSECTS(parcels.geom, luc.geom)
    AND
    ST_INTERSECTS(parcels.geom, wetlandyn.geom)
    AND
    CAST(parcels.camacls AS int) >=100
    AND
    CAST(parcels.camacls AS int) <200
    AND
    parcels.acres > 0
    AND
    ((ST_AREA(ST_INTERSECTION(ST_INTERSECTION(ST_INTERSECTION(parcels.geom,luc.geom),
    ST_INTERSECTION(parcels.geom,soils.geom)),ST_INTERSECTION(parcels.geom,wetlandyn.geom)))/ST_AREA(parcels.geom))*parcels.acres) 
    >= 0.00005

GROUP BY
    parcels.pin,
    soils.musym,
    parcels.camacls,
    luc.landcode,
    parcels.acres,
    parcels.geom,
    luc.geom,
    soils.geom,
    wetlandyn.geom,
    wetlandyn.inwetland

ORDER BY 
    parcels.pin, soils.musym, luc.landcode

I am not familiar with R-tree indexing at a source code level, but one thing I have noticed is that if you have a multipolygon that covers a wide area, compared to the polygons in your other tables, you might not get efficient use of the index, as all the polygons that intersect the bounding box of the multipolygon have to be checked for intersection -- this is my supposition. I have found that dumping multipolygons to individual polygons can dramatically speed up queries of this nature. — John Powell, Feb 20 '15 at 07:36

score 2 · Answer 1 · answered Feb 20 '15 at 01:33

2

The reason of why it takes so long for the query to run is because you are dealing with too many geometries at the same time. Perhaps you could create spatial indexes for your data, this usually greatly improves the processing times.

Here is some PostGIS documentation regarding indexes: http://postgis.net/docs/using_postgis_dbmanagement.html#idp7260032

Here is a quick guide of how to create spatial indexes in PostGIS: http://revenant.ca/www/postgis/workshop/indexing.html

answered Feb 20 '15 at 01:33

Ricardo Oliveira

283
3
14

Did the OP say there was no index? I somehow doubt that someone would write such a query and forget the index – John Powell Feb 20 '15 at 06:37
Fortunately, QGIS added the indexes when I imported them. I wouldn't have known how to do it myself until now. Thanks for the links. – Nathan Feb 20 '15 at 13:59

score 0 · Answer 2 · edited Apr 13 '17 at 12:34

0

I linked this up in Twitter and received this response from @pwramsey (https://twitter.com/pwramsey/status/568489615293845504)

"2 things: Acquiring ArcGIS-like speed in Postgis; and try exploding your wetlands into multiple singletons, not a few multis"

edited Apr 13 '17 at 12:34

Community

1

answered Feb 20 '15 at 13:53

Nathan

76
1
5

Well, I am not Paul Ramsey, but I am glad I said the same thing above. Split your multipolygons into single polygons, it makes the index work better, ie, avoid MBR on whole multpolygon. If you posted your explain, it is possible you will see a full table scan occuring, due to the multipolygon. – John Powell Feb 20 '15 at 16:22

Extracting data from four intersecting polygon layers in PostGIS?

2 Answers2