Follow up on this topic: Optimizing an intersection between a single massive multipolygon (WKT) and many features from PostGIS

I noticed I have to run an ANALYZE (not 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.

Source: https://www.2ndquadrant.com/en/blog/postgresql-vacuum-and-analyze-best-practice-tips/


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: (...)

Source: https://www.enterprisedb.com/blog/postgresql-query-optimization-performance-tuning-with-explain-analyze

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). https://registry.hub.docker.com/r/postgis/postgis/tags? landcover features splitted by the st_subdivded MultiPolygon
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)
  ST_Union( -- add ST_Union here to aggregate all splitted features
          WHEN ST_WITHIN(landcoverpolygons.geom, data.geom)
              THEN landcoverpolygons.geom
          ELSE ST_Intersection(data.geom, landcoverpolygons.geom)
  ) 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.

Using 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]     
    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"    
    WAGYU="0.5.0 (Internal)" TOPOLOGY
(1 row)

From postgis/postgis:13-3.1 docker image, available at:

  • 1
    What PostGIS version? This is likely fixed by GEOS 3.9 and up.
    – dr_jts
    Sep 29 at 16:16
  • Hmm, interesting. I'll give one of the available PG14-beta based image a try then, because I cannot figure out a PostGIS image based on PG13 which relies on a GEOS version > than 3.7.x (but I may have missed it). Sep 29 at 17:49
  • 1
    Results are simply outstanding with PG14! I feel like I've been propelled light years into the future! My profound respect to the developers. Really. Sep 29 at 20:35
  • 1
    Excellent news. GEOS 3.9 also implements an optimization for intersection with large geometries, so that might be helping here as well.
    – dr_jts
    Sep 30 at 18:46
  • 1
    One takeaway from this is that ST_Subdivide is probably not a great solution to improve intersection calculations. It's awkward to work with, and the new GEOS provides better performance anyway. It's still useful for improving performance of intersects queries, though.
    – dr_jts
    Sep 30 at 18:49

I've faced this error before, it likely has to do with polygon edges intersecting eachother and forming new polygons in the st_intersection function

  1. you can try and add st_makevalid to the st_intersection before it gets thrown through the st_union

    st_makevalid(ST_Intersection(data.geom, landcoverpolygons.geom))
  2. use st_collect and convert the collection to a polygon using st_buffer

           WHEN ST_WITHIN(landcoverpolygons.geom, data.geom) THEN landcoverpolygons.geom
           ELSE ST_Intersection(data.geom, landcoverpolygons.geom)
       AS intersect) ,0)::geometry(MultiPolygon,4326)
  3. if none of those work you can look into st_clusterdbscan (super powerful function but you will need to restructure a good portion of your code for this one)

  • 1
    It's more likely due to robustness issues in the old GEOS intersection code. That was dramatically improved in GEOS 3.9. The new implementation is much more robust, and also a lot faster for intersections between large polygons.
    – dr_jts
    Sep 30 at 18:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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