I noticed I have to run an
VACUUM ANALYZE!) on the table for the index to be effectively taken into account. Then,
ST_SubDivide really kicks in. This seems confirmed by:
The time when you must run ANALYZE manually is immediately after bulk loading data into the target table. A large number (even a few hundred) of new rows in an existing table will significantly skew its column data distribution. The new rows will cause any existing column statistics to be out-of-date. When the query optimizer uses such statistics, query performance can be really slow. In these cases, running the ANALYZE command immediately after a data load to completely rebuild the statistics is a better option than waiting for the autovacuum to kick in.
In the example above, the database had gone through a fair amount of activity, and the statistics were inaccurate. With an ANALYZE (not VACUUM ANALYZE or EXPLAIN ANALYZE, but just a plain ANALYZE), the statistics are fixed, and the query planner now chooses an Index Scan: (...)
But then, using ST_SubDivide leads to a first problem: if a landcover feature is covered by X parts of the subdivided MultiPolygon given as the input, it will results in X different records for that single landcover feature, where it should be one and only one result (Fig.1).
Fig.1 A single landcover feature is splitted by the ST_SubDivded MultiPolygon.
To fix that, I have to use ST_Union:
WITH data (geom) AS ( SELECT ST_Transform( ST_GeomFromText('MULTIPOLYGON (( <90'000 WGS84 points coordinates> )) ), 2056) ) SELECT id, name, ST_Union( -- add ST_Union here to aggregate all splitted features 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) -- group by each field in the select clause is mandatory ;
which, to me, will first probably annihilate the benefit of using
ST_SubDivide, but which also lead to this strange error:
ERROR: lwgeom_unaryunion_prec: GEOS Error: TopologyException: Input geom 1 is invalid: Self-intersection at or near point 8.546026915774 46.15793720531 at 8.546026915774 46.15793720531 SQL state: XX000
How could this be?
I'd like to know what may be the cause of this error, because it would be natural that if you clip some (multi)polygon, and aggregate the clipped parts, well... you should end up on the exact same original feature! But I face this strange GEOS bug here that I cannot understand. I have the feeling that it's related to some numerical rounding done under the hood by
ST_SubDivide, which makes this function maybe faster, but not reliable.
ST_Collect instead of
ST_Union is working, and as stated by the documentation of the latter:
ST_Collect may sometimes be used in place of ST_Union, if the result is not required to be non-overlapping. ST_Collect is usually faster than ST_Union because it performs no processing on the collected geometries.
But it results in "geometry collections", where I need single features.
Please note that applying
ST_Union on these geometry collections also raises the same error message as shown above.
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" GDAL="GDAL 2.4.0, released 2018/12/14" PROJ="Rel. 5.2.0, September 15th, 2018" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1" WAGYU="0.5.0 (Internal)" TOPOLOGY (1 row)
From postgis/postgis:13-3.1 docker image, available at: