Making slices through postgis raster tables

Question

I have an existing python script that i am trying to "translate" into postgres sql. The reason for this is that the server im using for this is not running python.

The objective is to slice through several postgis raster tabels intersecting a specific elevation value. Each raster is the bottom of a layer in a geological model. For a given elevation (e.g. 0 meter below sealevel) i want the cells where the elevation is between the bottom (current layer) and the top (the previous raster layer) displaying the layer id of the bottom layer.

My problem is that the sql is terrible and the time of execution is very long. I have 50 postgis raster tables and many elevations to loop through. In comparisson the python script, using numpy and rasterio, executed with half an hour for all elevations. This sql script uses the same time for a single elevation.

The figure below is my poor attempt to draw the objectiv. The left is displaying the input raster tables from sideview, as well as the elevation to slice them. The right displays the output raster from topview, displaying the layer id (e.g. layer 1 and 2) from the slice.

-- THIS IS ONLY RUNNING ONCE
DROP TABLE IF EXISTS mstfohm.table_name_out;

CREATE TABLE mstfohm.table_name_out AS 
WITH 
tmp AS 
(
    SELECT st_union(rast) AS rast
    FROM mstfohm.table_name_template
)
SELECT  
    ST_MapAlgebra(
            rast, 
            rast, 
            '[rast]^0'
    ) AS rast 
FROM tmp
;

-- THIS PART IS REPEATED A LOT OF TIMES FOR DIFFERENT ELEVATIONS AND RASTER TABLES
-- layer 1 
DROP TABLE IF EXISTS temp_raster;

CREATE TEMP TABLE temp_raster AS 
    (
        SELECT *
        FROM mstfohm.table_name_out
    )
;

DROP TABLE IF EXISTS mstfohm.table_name_out;

CREATE TABLE mstfohm.table_name_out AS 
WITH 
top AS 
(
    SELECT st_union(rast) AS rast
    FROM mstfohm.table_name_in1
),
bot AS 
(
    SELECT st_union(rast) AS rast
    FROM mstfohm.table_name_in2
),
topbot AS 
(
    SELECT 
        ST_MapAlgebra(
            top.rast, 
            bot.rast, 
            'CASE WHEN [rast1] >= 0 AND [rast2] < 0 THEN [rast1]^0 * 1 ELSE [rast2]^0 END'
            ) AS rast 
    FROM top, bot
)
SELECT 
    ST_MapAlgebra(
        a.rast, 
        b.rast, 
        'CASE WHEN [rast2] = 1 THEN [rast1] ELSE [rast2] END'
        ) AS rast
FROM topbot a, temp_raster b
;

-- layer 2 
DROP TABLE IF EXISTS temp_raster;

CREATE TEMP TABLE temp_raster AS 
    (
        SELECT *
        FROM mstfohm.table_name_out
    )
;

DROP TABLE IF EXISTS mstfohm.table_name_out;

CREATE TABLE mstfohm.table_name_out AS 
WITH 
top AS 
(
    SELECT st_union(rast) AS rast
    FROM mstfohm.table_name_in2
),
bot AS 
(
    SELECT st_union(rast) AS rast
    FROM mstfohm.table_name_in3
),
topbot AS 
(
    SELECT 
        ST_MapAlgebra(
            top.rast, 
            bot.rast, 
            'CASE WHEN [rast1] >= 0 AND [rast2] < 0 THEN [rast1]^0 * 2 ELSE [rast2]^0 END'
            ) AS rast 
    FROM top, bot
)
SELECT 
    ST_MapAlgebra(
        a.rast, 
        b.rast, 
        'CASE WHEN [rast2] = 1 THEN [rast1] ELSE [rast2] END'
        ) AS rast
FROM topbot a, temp_raster b
;

As a comparison the python script was in essence the code snippit below where the rasterio numpy array "array_slice" is updated through looping through all layers for each elevation.

array_slice = np.where(
    np.isnan(array_slice),
    np.where(
        (elevation_value >= array_bot) & (elevation_value <= array_top) (array_top != arr_bot),
        layer_id,
        array_slice 
    ),
    array_slice 
)