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

本文介绍了如何在PostGIS中使用自定义函数,通过数组操作和ST_MapAlgebra来实现每小时降雨量的栅格汇总,探讨了两种方法:一是将栅格数据转换为数组求和,二是循环调用ST_MapAlgebra进行逐帧累加。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

        近期有个需求要计算多幅栅格的总和, 比如: 气象降雨栅格每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(表达式版本)

数组函数和操作符

### 使用Python将栅格数据导入PostGIS数据库 要实现通过Python将栅格数据导入PostGIS数据库,可以采用种方式。以下是基于现有技术栈的一种解决方案。 #### 方法概述 一种常见的方式是使用GDAL/OGR工具结合Python来完成此操作。具体而言,可以通过调用`ogr2ogr`命令行工具或者直接利用GDAL库的功能来执行这一任务[^3]。此外,也可以借助SQL语句配合Python的`psycopg2`库手动插入栅格数据[^1]。 --- #### 方案一:使用 `ogr2ogr` 命令行工具(推荐) 这种方法简单高效,适合处理批量导入的任务。以下是一个完整的代码示例: ```python import os # 定义参数 host = "localhost" port = "5432" dbname = "your_db_name" user = "postgres" password = "your_password" input_raster_path = r"C:\path\to\your\raster_file.tif" # 构建 ogr2ogr 命令 command = f'ogr2ogr -f "PostgreSQL" PG:"host={host} port={port} dbname={dbname} user={user} password={password}" "{input_raster_path}"' # 执行命令 os.system(command) ``` 上述代码会将指定路径下的栅格文件导入到目标PostGIS数据库中。需要注意的是,输入的栅格文件应为支持的标准格式(如GeoTIFF),并且确保PostGIS扩展已启用于目标数据库[^4]。 --- #### 方案二:使用 GDAL 库直接读写 PostGIS 数据源 如果希望完全在Python环境中控制整个流程而不依赖外部命令行工具,则可考虑使用GDAL库。这种方式更加灵活但也稍显复杂。下面是具体的实现步骤: ```python from osgeo import gdal # 设置驱动程序名称 driver_name = "PostGISRaster" # 配置连接字符串 connection_string = ( f"PG:host={host} " f"port={port} " f"dbname='{dbname}' " f"user='{user}' " f"password='{password}'" ) # 注册所有驱动 gdal.AllRegister() # 打开原始栅格文件 source_dataset = gdal.Open(input_raster_path, gdal.GA_ReadOnly) if source_dataset is None: raise Exception(f"无法打开栅格文件 {input_raster_path}") # 创建新的PostGIS Raster数据集 target_dataset = gdal.GetDriverByName(driver_name).CreateCopy( connection_string, source_dataset, strict=0 ) del target_dataset # 关闭数据集以保存更改 ``` 在此过程中,需注意确认所使用的GDAL版本是否支持PostGIS Raster驱动以及相应的配置选项[^3]。 --- #### 方案三:通过 psycopg2 插入栅格数据 另一种可能的方法涉及直接构建SQL查询并通过`psycopg2`将其发送至服务器端执行。这通常适用于单张图像的小规模场景下。例如: ```python import psycopg2 from rasterio.transform import Affine from rasterio.crs import CRS # 连接到数据库 conn = psycopg2.connect( host="localhost", port="5432", user="postgres", password="your_password", dbname="your_db_name" ) cur = conn.cursor() # 准备栅格数据及其元信息 with open(input_raster_path, 'rb') as f: byte_data = f.read() transform_matrix = (Affine(1.0, 0.0, 0.0, 0.0, -1.0, 0.0)).to_gdal()[:6] srid_value = int(CRS.from_epsg(4326).to_authority()[1]) # 插入记录 query = """ INSERT INTO your_table_name (rast) VALUES (ST_SetSRID(ST_MakeEmptyRaster(%s, %s, %s, %s, %s, %s, %s, %s), %s)); """ params = (*transform_matrix, srid_value) cur.execute(query, params) # 提交事务并关闭资源 conn.commit() cur.close() conn.close() ``` 请注意替换模板中的占位符以匹配实际需求,并验证表结构定义是否兼容预期的操作模式[^1]。 --- ### 总结 以上介绍了三种不同的策略用于解决如何运用Python向PostGIS数据库上传栅格资料的问题。每种方案各有优劣,在选用前建议综合考量项目背景、环境约束和技术偏好等因素做出决定。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值