Using ST_Union to aggregate polygon parts resulting of ST_SubDivide leads to a GEOS TopologyException

Question

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/

And:

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?
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.

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.

Versioning:

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:
https://registry.hub.docker.com/r/postgis/postgis/tags

Answer 1

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

   ST_Buffer(ST_Collect(
      CASE
          WHEN ST_WITHIN(landcoverpolygons.geom, data.geom) THEN landcoverpolygons.geom
          ELSE ST_Intersection(data.geom, landcoverpolygons.geom)
      END
      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)