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

Question

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

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

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

Answer 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;

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