Optimizing an intersection between a single massive multipolygon (WKT) and many features from PostGIS

Context

I have a massive MultiPolygon (~90'000 vertices, super complex, spreading over 60km², with lots of holes and parts) resulting of some Python computation. This MultiPolygon is coming as a WKT string. When I copy paste the query with that MultiPolygon into pgAdmin4, it makes the browser freeze for ~1minute. And if I dump that query into a file, it takes some time, and the file is ~20MB.

On the other hand, I have some other polygons (~86'000 features). They are already stored in a PostGIS database. They correspond to some specific regions of interest over a country. They are randomly spread out all over the country. They do not form a continuous surface. And they are quite small (to some square meters to max 1-2km²).

Needs

I need to get all features which intersect the large MultiPolygon. I am especially interested in the centroids of the polygons, after being clipped by the fat MultiPolygon (the ones outside being naturally no more needed).

Issue and SQL Code

Here's the query, it's embedded as a sub-sub `SELECT` statement (I found it was faster to "split" the code like this). I identified the part shown here as the bottleneck. It's taking ~1h on a mainstream 2016 laptop.

``````SELECT  fid, id, name, stuff, pseudo_d3d
FROM (
SELECT
*,
round(ST_3ddistance(
ST_SetSRID(st_makepoint(2600000.0,1200000.0,500.0), 2056),
ST_SetSRID(
ST_MakePoint(ST_X(centroid),ST_Y(centroid),ST_Z(intersect)/0.9989
), 2056)
)
) AS pseudo_d3d
FROM (
SELECT
ST_Intersection( -- compute the real intersection
landcoverpolygons.geom,
SELECT ST_Transform(
ST_GeomFromText('MULTIPOLYGON (( <90'000 WGS84 points coordinates> ))
), 2056)
) AS intersect
FROM landcoverpolygons
-- polygons_categories is holding some values list
INNER JOIN polygons_categories
ON (landcoverpolygons.landcovertype = polygons_categories.landcovertype)
WHERE ST_Intersects( -- compute the intersect to filter out prior to the intersection
ST_Transform(
landcoverpolygons.geom, 2056
),
SELECT ST_Transform(
ST_GeomFromText('MULTIPOLYGON (( <90'000 WGS84 points coordinates> ))
), 2056)
)
AND landcoverpolygons.landcovertype IN (
SELECT landcovertype FROM polygons_categories WHERE forest=True
)
AND landcoverpolygons.name_type <> 'foreign'
) AS foo LATERAL ST_SetSRID(ST_Centroid(foo.intersect), 2056) AS centroid
) AS bar;
``````

What I've tried so far

A first idea was to test for the `ST_Intersects` prior to the computation of `ST_Intersection`. This improved things consequently (between ~50 and 400% depending on the input MultiPolygon).

In addition, as the following part repeats two times...:

``````SELECT ST_Transform(
ST_GeomFromText('MULTIPOLYGON (( <90'000 WGS84 points coordinates> ))
), 2056)
``````

...I would wrote this as follow using a `WITH` statement if that query wasn't enclosed in two other `SELECT` statements:

``````WITH data (foo) AS (
SELECT ST_Transform(
ST_GeomFromText('MULTIPOLYGON (( <90'000 WGS84 points coordinates> ))
), 2056)
)
SELECT
ST_Intersection( -- compute the real intersection
landcoverpolygons.geom,
data.foo
) AS intersect
FROM landcoverpolygons
INNER JOIN polygons_categories
ON (landcoverpolygons.landcovertype = polygons_categories.landcovertype)
WHERE ST_Intersects( -- compute the intersect to filter out prior to the intersection
ST_Transform(
landcoverpolygons.geom, 2056
),
data.foo
)
AND landcoverpolygons.landcovertype IN (
SELECT landcovertype FROM polygons_categories WHERE forest=True
)
AND landcoverpolygons.name_type <> 'foreign';
``````

This seems to go "a little" faster (~<5%).

Please, also notice that I already have indexes on `landcoverpolygons.geom`, `landcoverpolygons.landcovertype` and `polygons_categories.landcovertype`. I also have indexes on the transformation shown here.

So, indexes on `landcoverpolygons` are as follow:

``````CREATE INDEX landcoverpolygons_wgs84_gix
ON public.landcoverpolygons
USING gist (geom);

CREATE INDEX landcoverpolygons_2056_gix
ON public.landcoverpolygons
USING gist (ST_Transform(geom, 2056));
``````

But it's not at all improving the computation time. It's even the opposite and made it rise between +50% and +100%! Surprisingly. Or not.... Maybe because there are features all over the massive MultiPolygon, so I would intuitively say that the index is of absolutely no help here (which seems confirmed here: https://stackoverflow.com/a/5203827/6630397).

Question

How could I lower the computation time of `ST_Intersect()` which seems to be the most time consuming operation here (see edits below)?

EDIT:

• As suggested by geozelot, using `ST_SubDivide()` improved the computation time between 0-6% max.: 56 minutes instead of 1h in the best run I've executed. In addition, I faced this error when using `ST_SubDivide()` in the `WHERE` statement, so I was only able to use it when first using the MultiPolygon (in the first version of the code, without the CTE, which is the one I use when actually running the full query):
``````set-returning functions are not allowed in WHERE
``````

In addition, the problem with ST_SubDivide() is that if a landcover feature is covered by X parts of the subdivided MultiPolygon, it will results in X different records for that exact same landcover feature, but I naturally only need a single one (Fig.1).

Fig.1 As a result, a single landcover features is splitted by the st_subdivded MultiPolygon.

• As suggested by ziggy, I also tried to replace `data.foo` in the ST_Intersect block (but not the ST_Intersection of course!) by `ST_Envelope(data.foo)` or even `ST_ConvexHull(data.foo)` but this also doesn't improve the computation time.

From these observations, I think I can pretty much say that it's the `ST_Intersections()` itself which costs the most time.

EDIT 2021-09-26: I also discovered the 'EXPLAIN ANALYZE' button on pgAdmin4 (but I don't understand well how to interpret the results, this is quite obscure to me).

I think this, from the 'Analysis' tab in the bottom pane showing the results may be relevant (this is another feature I'm actually working on, which only takes some minutes instead of a full hour):

Fig.2. Analysis tab in the results pane of pgAdmin4 after EXPLAIN ANALYZE was run with all options checked.

along with this table, from the 'Statistics' pane:

Fig.3 "Statistics per Node Type" from the 'Statistics tab in the pgAdmin4 results pane.

But this surprises me as the shown condition is not mentioning the intersection. And I have no idea what is this `_1` suffix doing (second row of Fig.2).

Versioning:

PostgreSQL: "PostgreSQL 13.3 (Debian 13.3-1.pgdg100+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 8.3.0-6) 8.3.0, 64-bit"
PostGIS: "3.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1"

• Hey this is a good question(s), but it's 3 questions in one. Can you focus on one issue at a time (maybe ask multiple questions)? Also, can you share the rest of your code (the ... parts)? You say 'I need to get the smaller polygons that intersect the large polygon'. You don't need st_intersection for that. Do you also need to get the actual intersection between the polygons? One very obvious problem I see is that you're transforming all of your geometries. Unless you have a transformation index on both tables, the spatial index won't be used on either table. Sep 24, 2021 at 18:27
• `ST_SubDivide` your MultiPolygon (e.g. in a CTE) Sep 24, 2021 at 19:11
• i suggest forcing a bounding box intersect before any of the hard geom spatial functions Sep 25, 2021 at 2:48
• i also suggest using a case statement and st_within if any of the polygons fall within the foo data Sep 25, 2021 at 2:51
• Please bear in mind that throwing indexes at everything is not helpful. In theory, an index on the main decision column on the larger of two related tables is the most performant combination - while it would also work the other way around. In every case, indexes are considered for only one of two relations, i.e. one per joined relation. If your base table has an index that can be used, and is the running table, no need for another index elsewhere. Sep 25, 2021 at 9:30

okay a few things as I mentioned in the comments

1. on larger intersection tables its always faster to force a bounding box intersect to weed out non intersections using &&

2. you want to avoid st_transform in the where clauses, it can slow it down (so the query I post is going to disregard your st_transform on the fly for the large table because you shouldnt be doing that if you could avoid it

3. Sometimes putting the attribute where clauses before the spatial ones could speed up the query

4. I added a case statement to the top portion of the st_intersection. you want to avoid using this function if possible because its so expensive. So if you know any of these polygons will fall into each other completely then just grab that entire geometry instead of running an unnecessary intersection

5. im not touching your main outer query because im not entirely sure why you have a lateral there and that doesnt seem to be the source of the slowness

`````` WITH data (foo) AS (
SELECT ST_Transform(
ST_GeomFromText('MULTIPOLYGON (( <90'000 WGS84 points coordinates> ))
), 2056)
)

SELECT  fid, id, name, stuff, pseudo_d3d
FROM (
SELECT
*,
round(ST_3ddistance(
ST_SetSRID(st_makepoint(2600000.0,1200000.0,500.0), 2056),
ST_SetSRID(
ST_MakePoint(ST_X(centroid),ST_Y(centroid),ST_Z(intersect)/0.9989
), 2056)
)
) AS pseudo_d3d
FROM (
SELECT
case -- need to asses if the landcoverpolygons fall within data or vice versa, if vice versa then switch the inputs
WHEN ST_WITHIN(landcoverpolygons.geom,data.foo) then landcoverpolygons.geom
else  ST_Intersection(landcoverpolygons.geom,  data.foo)
end AS intersect
FROM landcoverpolygons
JOIN polygons_categories
ON (landcoverpolygons.landcovertype = polygons_categories.landcovertype)
WHERE  landcoverpolygons.landcovertype IN (
SELECT landcovertype FROM polygons_categories WHERE forest=True
)
AND landcoverpolygons.name_type <> 'foreign'
AND landcoverpolygons.geom && data.foo
AND ST_Intersects(landcoverpolygons.geom, data.foo)
) AS foo LATERAL ST_SetSRID(ST_Centroid(foo.intersect), 2056) AS centroid
) AS bar;
``````
• @swiss_knight you don't have tmpdata in your example? Sep 25, 2021 at 23:40
• Thanks, Will try that! But a small detail that I can see right now is that we also need to call `data` in the latest `FROM` clause, so: `FROM landcoverpolygons` must be `FROM landcoverpolygons, data` otherwise the `ERROR: missing FROM-clause entry for table "data"` will be raised. Sep 26, 2021 at 6:38
• it was an error due to copy/pasting an old test I did... Previous comment fixed to align to what's in the thread. Sep 26, 2021 at 7:23

Point set intersections are very costly, and you want to absolutely make sure you limit your queries

• to the least possible function calls
• on the least possible amount of vertices between both participating geometries

For this you want to

• first `ST_SubDivide` your MultiPolygon, resulting in some 300 times less vertices per operation
• then explicitly `JOIN` on `ST_Intersects` to make use of the now more effective index on `landcoverpolygons.geom`, limiting operations to only the relevant subset
• while using a `CASE` decision to avoid any unnecessary call to `ST_Intersection`
• and finally get your `ST_Centroid` on the `ST_Collect` of those individual intersections that you mention above

Note that

• it doesn't make sense to transform your geometries at all - especially if `landcoverpolygons.geom` is actually referenced in `EPSG:4326`

• there is no benefit in any of your spatial indexes if more than a few percent of rows would match, and even if the planer chose to use it over the attribute indexes on your relations, it would be highly ineffective, producing a high number of false hits

• Multi geometries are always bad for spatial relationship analysis, and even more so are geometries with excessive amounts of vertices

• adding indexes by guess is likely having the opposite effect - get rid of all excessive indexes, run `VACUUM ANALYZE` on all participating tables and retry; for this particular query only the indexes

• `ON landcoverpolygons USING GIST (geom)`
• `ON landcovertype (forest)`

are of interest

Try:

``````WITH
mp_sd AS (
SELECT  ST_Subdivide(<MP>) AS geom
)
SELECT  id,
ST_Centroid(its_geom) AS geom
FROM    (
SELECT  lcp.id,
ST_Collect(
CASE ST_Within(lcp.geom, mp_sd.geom)
WHEN TRUE THEN lcp.geom
ELSE ST_Intersection(lcp.geom, mp_sd.geom)
END
) AS its_geom
FROM    landcoverpolygons AS lcp
JOIN    mp_sd
ON    lcp.geom && mp_sd.geom AND ST_Intersects(lcp.geom, mp_sd.geom)
WHERE   lcp.name_type <> 'foreign'
AND   lcp.landcovertype IN (
SELECT  landcovertype
FROM    polygons_categories
WHERE   forest
)
GROUP BY
1
) q
;
``````

Always post the full `EXPLAIN ANALYZE` output here - some of us know how to make sense of it and may be able to help you even more!

• good call with the st_subdivide Sep 28, 2021 at 14:33

Previous answers are of course precious, but what I discovered by testing a PG14 based image is... amazing! I cannot believe my eyes after some early testing (maybe this is simply related to a better default configuration, bfut if not, improvements are really outstanding)!

I took the exact same query that was problematic to me on PG13+PostGIS3.1:

• with spatial indexes:
``````Successfully run. Total query runtime: 3 min 4 secs.
207 rows affected.
``````
• without spatial indexes:
``````Successfully run. Total query runtime: 1 min 2 secs.
207 rows affected.
``````
• with/without spatial indexes but with the `ST_SubDivide()` + `ST_Union()` trick:
buggy

Same query copied/pasted it into PG14+PostGIS3.1:

• with spatial indexes:
``````Successfully run. Total query runtime: 10 secs 831 msec.
207 rows affected.
``````

(yes, 10 seconds!)

• without spatial indexes:
``````Successfully run. Total query runtime: 9 secs 979 msec.
207 rows affected.
``````
• with spatial indexes but with the `ST_SubDivide()` + `ST_Union()` trick:
``````Successfully run. Total query runtime: 1 min.
207 rows affected.
``````
• without spatial indexes but with the `ST_SubDivide()` + `ST_Union()` trick:
``````Successfully run. Total query runtime: 42 secs 609 msec.
207 rows affected.
``````

/Same machine, freshly build container in both situations.

The request was:

``````WITH data (geom) AS (
SELECT
ST_SubDivide --<<
(
ST_Transform(
ST_GeomFromText('MULTIPOLYGON (( <90'000 WGS84 points coordinates> ))
), 2056)
)
)
SELECT
id,
name,
ST_Union --<<
(
CASE
WHEN ST_WITHIN(landcoverpolygons.geom, data.geom)
THEN landcoverpolygons.geom
ELSE ST_Intersection(data.geom, landcoverpolygons.geom)
END
)
AS intersect
FROM data, landcoverpolygons
INNER JOIN polygons_categories
ON (polygons_categories.landcovertype = landcoverpolygons.landcovertype)
WHERE landcoverpolygons.landcovertype
IN (SELECT landcovertype FROM polygons_categories WHERE forest=True)
AND landcoverpolygons.name_type <> 'foreign'
AND landcoverpolygons.geom && data.geom
AND ST_Intersects(landcoverpolygons.geom, data.geom)

GROUP BY (name, id) --<<
;
``````

The 3 lines with `--<<` are the ones that I added when I said: "with the `ST_SubDivide()` + `ST_Union()`".

Versioning:

PG13: from postgis/postgis:13-3.1

``````SELECT version();
----------------------------------------
PostgreSQL 13.3 (Debian 13.3-1.pgdg100+1)
on x86_64-pc-linux-gnu,
compiled by gcc (Debian 8.3.0-6) 8.3.0, 64-bit
(1 row)

SELECT postgis_full_version();
----------------------------------------
POSTGIS="3.1.2 cbe925d" [EXTENSION]
PGSQL="130"
GEOS="3.7.1-CAPI-1.11.1 27a5e771"
PROJ="Rel. 5.2.0, September 15th, 2018"
GDAL="GDAL 2.4.0, released 2018/12/14"
LIBXML="2.9.4"
LIBJSON="0.12.1"
LIBPROTOBUF="1.3.1"
WAGYU="0.5.0 (Internal)" TOPOLOGY
(1 row)
``````

PG14: from postgis/postgis:14rc1-3.1

``````SELECT version();
----------------------------------------
PostgreSQL 13.3 (Debian 13.3-1.pgdg100+1)
on x86_64-pc-linux-gnu,
compiled by gcc (Debian 8.3.0-6) 8.3.0, 64-bit
(1 row)

SELECT postgis_full_version();
----------------------------------------
POSTGIS="3.1.3 008d2db" [EXTENSION]
PGSQL="140"
GEOS="3.9.0-CAPI-1.16.2"
PROJ="7.2.1"
GDAL="GDAL 3.2.2, released 2021/03/05"
LIBXML="2.9.10"
LIBJSON="0.15"
LIBPROTOBUF="1.3.3"
WAGYU="0.5.0 (Internal)" TOPOLOGY
RASTER
(1 row)
``````

I'll let you judge for yourself.

On Thursday, 11 April 2019, a SpaceX Falcon Heavy rocket launched the Arabsat-6A satellite
from Launch Complex 39A at NASA’s Kennedy Space Center in Florida.

Source

• That's great. Can you post the query that was fastest? Sep 30, 2021 at 18:42
• Does the faster query use ST_Intersection on the original large polygon (as opposed to the ST_Subdivide approach)? Sep 30, 2021 at 18:50
• As I noted on the other issue, GEOS 3.9 implements an optimization for intersection with large geometries. That is probably what is improving performance here. Sep 30, 2021 at 19:05
• What is the `select postgis_geos_version()` in those images? Might be the main difference. Sep 30, 2021 at 20:36
• @dr_jts query added. Oct 2, 2021 at 8:13