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

Question

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"

Answer 1

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;
   

Answer 2

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!

Answer 3

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