ST_3DDifference in PostGIS taking several minutes to run on two POLYHEDRALSRUFACES

Question

I would like to calculate the 3D volume of a dam in a valley using PostGIS.

dam_polysurf is created as a POLYHEDRALSURFACE by extruding a vertical triangular polygon along the axis of the dam.

valley_polysurf is created as follows:

SELECT ST_Delaunaytriangles (ST_Collect ( ARRAY ( SELECT ( (ST_Dumppoints ( ST_Intersection (mcs.geom, st_buffer(da.geom, 100)) ) ).geom ) from "10m_contours_simp2" mcs, dam_axis da ) ), 0.0, 2 ) geom INTO test_tin2

SELECT ST_Extrude(tt.geom, 0, 0, -100) geom INTO valley_polysurf FROM test_tin2 tt

I then use ST_3DDifference to calculate the volume of the dam above the valley polysurface. The problem I have is that this takes a long time:

EXPLAIN ANALYZE SELECT ST_3ddifference(dp.geom, vp.geom) FROM dam_polysurf dp, valley_polysurf vp

Nested Loop (cost=0.00..485570.60 rows=1849600 width=32) (actual time=579640.077..579640.083 rows=1 loops=1)

-> Seq Scan on dam_polysurf dp (cost=0.00..23.60 rows=1360 width=32) (actual time=0.025..0.025 rows=1 loops=1)

-> Materialize (cost=0.00..30.40 rows=1360 width=32) (actual time=0.030..0.035 rows=1 loops=1)

-> Seq Scan on valley_polysurf vp (cost=0.00..23.60 rows=1360 width=32) (actual time=0.025..0.028 rows=1 loops=1) Planning Time: 0.560 ms Execution Time: 579640.261 ms

I don't know how to get the number of faces directly from a POLYHEDRALSURFACE. However using:

SELECT (st_dump(geom)).geom FROM valley_polysurf = 767 POLYGONZ

SELECT (st_dump(geom)).geom FROM dam_polysurf = 5 POLYGONZ

The planner thinks there are a lot of row involved so I have tried adding a gist_geometry_ops_nd index to both tables. However as they only have one row each this does not improve the actual query time.

I have tried clipping the base to the valley_polysurf to a constant z-value which reduces the number of POLYGONZ to 544, reducing the ST_3ddifference time but the intermediate process of clipping (also using ST_3ddifference) takes more time than is saved.

I assume that the algorithm is doing some kind of CROSS JOIN triangle-wise test to see if they intersect to form a line like this, and then creating a new set of surfaces for the diffrence. I have checked the SFCGAL website but it seems you have to make the documentation yourself (I don't have Doxygen installed) so I cannot confirm this.

Using:

POSTGIS="3.5.2 3.5.2" [EXTENSION] PGSQL="150" GEOS="3.13.0-CAPI-1.19.0" SFCGAL="SFCGAL 1.5.2, CGAL 5.6.1, BOOST 1.84.0" PROJ="8.2.1

Any thoughts on how to improve the execution time?