如何平均单个 PostGIS 栅格表中的所有波段？

Posted 2023-03-24

【中文标题】如何平均单个 PostGIS 栅格表中的所有波段？

【英文标题】：How can I average all bands in a single PostGIS raster table?

【发布时间】：2021-05-10 04:50:29

【问题描述】：

我在启用了 postgis 的 postgres 数据库中的单个栅格表中有 1 年的每日气候数据。它有 365 个乐队（每天一个）。如何平均所有波段以获得每个像素的单个年平均值？我试过ST_Union，但它返回了所有波段，或者我没有正确使用它：

select rid, st_union(rast, 'MEAN')
from climate_table
group by rid;

我想出了一个使用 ST_DumpAsPolygons 的解决方法，但它非常慢。任何建议表示赞赏。 （另外，我不敢相信“乐队”还不是一个标签，而且我没有足够的声誉来创建它！）

【问题讨论】：

ST_Union 如果您每天有一个光栅，每个光栅有 1 个波段，则应该可以工作。 （每天一分贝行）。您可能可以将 generate_series(1,365) 用于波段索引并使用 ST_Band(original_raster,index) 加入该索引，然后使用 ST_Union(rast,'MEAN') 生成的栅格集

@clamp，我有一个大栅格（约 500,000 个单元格），每个波段都是单独的一天。如果我有“每天一个光栅，每个光栅有 1 个波段（每天一个 db 行）”，我就不会有这个问题。此外，“每天一个光栅，每个 1 个波段”和“每天一个 db 行”听起来像是 2 个非常不同的东西

该光栅文件是否可供公众访问？

ST_Union 采用一组栅格 - 而不是一组波段，我提到了可用于在查询中转换栅格的函数。您想要一个产生 365 行且每行一个栅格（即每天）的 with 子句或子选择。这可以通过 ST_Union 聚合

可以从这里获得作为 NetCDF 的栅格：gdo-dcp.ucllnl.org/pub/dcp/archive/cmip5/loca/LOCA_2016-04-02/… 我实际上需要处理大约 1,000 个栅格，因此我在寻找一种更快的方法。感谢您查看这个。我会试试你的测试用例，但我开始认为它不会比使用 ST_DumpAsPolygons 或 ST_PixelsAsPoints 快。我会及时通知你。

【参考方案1】：

根据the documentation，函数ST_Union() 采用一组栅格。您有一个包含许多波段的栅格。 您可以使用ST_Union，但您必须将数据转换为一组栅格，其中每个栅格首先包含一个波段。

这是一个测试用例：

-- Create a table and add a row with an empty raster:

CREATE TABLE r_src(rid integer primary key, rast raster);
INSERT INTO r_src values(1,ST_MakeEmptyRaster(4,4,0,0,1));

-- Add two band layers to this raster one with value 10 for each point,
-- the other with value 4 for each point.

UPDATE r_src
    SET rast = ST_AddBand(rast,'4BUI'::text,10)
WHERE rid = 1;
UPDATE r_src
    SET rast = ST_AddBand(rast,'4BUI'::text,4)
WHERE rid = 1;

-- check what we have

SELECT  (rmd).width, (rmd).height, (rmd).numbands
FROM (SELECT ST_MetaData(rast) As rmd
    FROM r_src WHERE rid = 1) AS foo;
-- 4| 4| 2  --ok

-- output rasterdata in hex for both bands:

SELECT ST_AsHexWKB(rast) As rastbin FROM r_src WHERE rid=1;


-- output rasterdata for the 2nd band (note the '...040404' point values 
-- compared to '...0A0A0A' if you select the band at idx 1

SELECT ST_AsHexWKB(ST_Band(rast,2))FROM r_src WHERE rid = 1;


-- now step by step...
-- src is the single rast in our source table,
-- in extracted this is joined with a generated series
-- for each value in the series a new raster is returned 
-- containing the band at that index:

WITH src       AS (SELECT rast FROM r_src WHERE rid = 1),
     extracted AS (SELECT ST_Band(src.rast,ser.idx) 
                   FROM src
                   JOIN generate_series(1,2) AS ser(idx ) ON true)
SELECT * FROM extracted;

-- outputs 2 rows with one raster each, and each raster containing one band

-- now aggreagating these:

WITH src       AS (SELECT rast FROM r_src WHERE rid = 1),
     extracted AS (SELECT ST_Band(src.rast,ser.idx) AS r_day
                   FROM src
                   JOIN generate_series(1,2) AS ser(idx ) ON true)
SELECT ST_Union(r_day,1,'MEAN') FROM extracted;

-- outputs one row with the resulting raster holding one band 
-- that band contains the 'MEAN' point values '...070707'

我使用您 linked 的数据集检查了这一点，并使用 raster2pgsql 将其导入到 postgres 中。由于它的大小，我决定将它分成 100x100 的图块，并得到一个 45 行的表。现在要转换整个栅格，我使用了以下语句：

CREATE TABLE averages AS
WITH extracted AS (SELECT src.rid , ST_Band(src.rast,ser.idx) AS r_day
                   FROM rastertest src
                   JOIN generate_series(1,365) AS ser(idx ) ON true)
SELECT rid, ST_Union(r_day,1,'MEAN') AS rast
FROM extracted
GROUP BY rid;

这将创建一个仅包含平均值的新表平均值。

【讨论】：

非常感谢@clamp 看到这一点。对不起我之前的困惑。你解释得很好，你的例子很完美。我为我的特定问题发现的另一种可能的方法是在加载到 db 之前在 netcdf 文件上使用 gdal_calc.py。我还发现使用大图块（即 480x245，或只有 4 行）可以非常快速地加载这些栅格，并且 ST_Union 函数花费相同的时间 - 在加载/处理数百个这些时很有用。

感谢您的反馈。创建平均值的语句在此处完成需要 3 分钟（对于整个栅格）。我还认为在导入之前这样做可能会更快......