通过 postgis 栅格表制作切片

Question

我有一个现有的 python 脚本，我正在尝试将其“翻译”为 postgres sql。 原因是我使用的服务器没有运行 python。

目标是对与特定高程值相交的多个 postgis 栅格表进行切片。每个栅格都是地质模型中图层的底部。对于给定的海拔（例如，海平面以下 0 米），我希望海拔位于底部（当前图层）和顶部（前一个栅格图层）之间的单元格显示底层的图层 ID。

我的问题是sql很糟糕，执行时间很长。我有 50 个光栅文件和许多高程要循环。相比之下，Python 脚本使用 numpy 和 rasterio，对于所有海拔高度执行时间为半小时。此 SQL 脚本对单次海拔使用相同的时间。

下图是我绘制目标的拙劣尝试。左侧显示侧视图中的输入栅格表，以及对它们进行切片的高程。右侧显示俯视图的输出栅格，显示切片的图层 ID（例如图层 1 和图层 2）。

-- 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, 
            'CASE WHEN [rast]<0 THEN [rast]^0 ELSE [rast]^0 END'
    ) 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
;