PostGis 多栅格数据求和函数实现

        近期有个需求要计算多幅栅格的总和, 比如: 气象降雨栅格每5分钟更新一幅,每个小时有12幅栅格,业务要求统计每个小时的降雨量, 这就需要把这12幅栅格数据进行求和汇总为一幅新栅格以便于后续处理.不知道PostGis是否已有相关的函数, 如果有的话还请不吝赐教.由于未找到可用函数,隧而使用自定义函数处理.

        两个思路:

        ①提取栅格数据为二维数组,对数据进行求和,利用ST_SetValues函数更新栅格数据

        ②编写函数, 循环调用ST_MapAlgebra函数进行汇总求和, 直接得到最终栅格

        函数如下:

        数组计算方式:

主要函数:

1, 二维数组转一维数组(在用)

CREATE OR REPLACE FUNCTION unnest_2d_1d(ANYARRAY, OUT a ANYARRAY)
  RETURNS SETOF ANYARRAY AS
$func$
BEGIN
   FOREACH a SLICE 1 IN ARRAY $1 LOOP
      RETURN NEXT;
   END LOOP;
END
$func$ LANGUAGE plpgsql IMMUTABLE STRICT;

2, 一维数组转二维数组(未用) 

CREATE OR REPLACE FUNCTION "public"."array_agg_2d"("arrin" anyarray, "xsize" int4)
  RETURNS "pg_catalog"."anyarray" AS $BODY$
DECLARE
    colsn integer := 0;
    r RECORD;
    arr numeric[] := array[]::numeric[];
    arrout numeric[][] := array[]::numeric[][];
BEGIN
    FOR r IN (SELECT UNNEST(arrin) AS val) LOOP
      colsn = colsn + 1;
        select array_append(arr, r.val) into arr;
          IF colsn = xsize THEN
              select array_agg(arr) into arrout;
                arr := array[]::numeric[];
            ELSEIF colsn % xsize = 0 THEN
              select array_cat(arrout, arr) into arrout;
          arr := array[]::numeric[];
                colsn :=0;
            END IF;    
  END LOOP;
    RETURN arrout;
END
$BODY$
  LANGUAGE plpgsql IMMUTABLE STRICT
  COST 100;

示例代码: 

with arr2 as (
    SELECT array_agg(val)::numeric[][] as arr
    from (
        SELECT array_agg(val) val
        from (
            select rn, sum(arr[idx]) val
            from (
                select case row_number() over()%height when 0 then 9999 else row_number() over()%height end as rn, arr
                from (
                    select st_height(rast) as height, unnest_2d_1d(st_dumpvalues(rast,1)) as arr from g_rast
                ) dim1
            ) t, generate_subscripts(t.arr, 1) AS idx
            group by rn, idx
            order by rn, idx
        ) rs
        group by rn
        order by rn
    ) dim2
)
select st_setvalues(st_addband(ST_MakeEmptyRaster(rast), '32BF'::text, 0), 1, 1, 1, arr::numeric[][]) 
from g_rast r, arr2 where r.id = 1;

栅格计算方式

主要函数: 

1, 双栅格求和
CREATE OR REPLACE function sum_raster_state(raster, raster)
returns raster
language sql
as $f$

SELECT
CASE $1
WHEN NULL THEN
      $2
ELSE
    ST_MapAlgebra(
        $1, 
        $2, 
        '([rast1] + [rast2])', 
        NULL, 
        'UNION', 
        '[rast2]', 
        '[rast1]', 
        NULL)
END;
$f$;

2, 单栅格处理

CREATE OR REPLACE FUNCTION sum_raster_final(raster)
returns raster
language sql
as $f$
SELECT 
   $1
$f$;

3, 定义聚合函数

create aggregate sum_raster(raster) (
    SFUNC = sum_raster_state,
    STYPE = raster,
    FINALFUNC = sum_raster_final
);   

代码示例:

select sum_raster(rast) from g_rast

主要参考连接:

PostGIS栅格求和（地图代数） 

ST_MapAlgebra(表达式版本)

数组函数和操作符