实战课堂 | 科技减灾如何使用Geo-SQL实现洪灾承灾体损失评估？

2021-09-17 01:17:09

????1. 简介

洪涝灾害是我国目前面临的最主要的自然灾害，利用洪涝灾害承灾体损失综合评估模型，对灾害损失率和损失价值分布进行科学地计算，对于指导洪涝救灾、建立灾害预警机制、加强洪涝灾害成灾机制的研究，建立和完善更科学、更准确的洪灾损失评估预测体系具有重要的意义。

????1.1 模型简介

该模型利用洪水灾害淹没水深分布，结合承灾体类型、承灾体价值及脆弱性数据，计算灾害损失率和损失价值分布。

1.1.1 数据输入

输入的数据包括：

洪水淹没水深分布
土地利用类型分布
土地利用价值分布
土地利用淹没水深-损失率对照表

1.1.2 数据计算

1，承灾体损失率计算公式：
DR=f(H)

其中DR为每类承灾体的损失率；H为致灾因子强度；f为致灾因子和损失率映射关系。

2，承灾体损失价值计算公式：
Loss = DR * E

其中Loss为损失价值，DR为灾害损失率，E为承灾体价值。

1.1.3 数据输出

洪涝灾害损失率分布
洪涝灾害损失价值分布

????1.2 Ganos

Ganos是阿里云自研的时空数据库引擎，包含了几何，栅格，点云，几何网络和轨迹模型5大数据模型，支持RDS POSTGRESQL 和 POLARDB-P*品。
其中针对时空栅格数据类型，Ganos提供了超大规格栅格数据存储、计算能力，单幅栅格数据理论上无容量限制，具备全球一张图的管理能力，使得传统GIS中复杂的栅格分析操作可以使用Geo-SQL轻松实现，并具备了与几何数据类型统一分析能力。

具体函数见 ????https://help.aliyun.com/document_detail/107567.html
地图代数是栅格分析与GIS建模中常用技术方法。Ganos为栅格图层计算操作提供了栅格代数表达式语言和一组栅格代数函数，称为ACL(Algebra Computing Language)。

ACL包括通用算术运算符，逻辑运算符，位运算，关系运算符以及一组统计函数，并允许它们*组合使用，实现更为复杂的运算操作。Ganos 栅格借助于ACL强大的计算表达式，支持基于一个或多个栅格对象像元值的条件查询，数学建模，分类操作以及生成新的结果栅格对象。

PL/pgSQL与ACL结合使用，则提供了更为强大的易用的栅格分析工具。PL/pgSQL提供变量和常量的声明，通用数学表达式，基本函数，逻辑判断和流程控制，ACL为栅格计算提供像元代数计算的表达方式。用户可以轻松地结合两种的优势进行时空栅格的分析与建模，如对全球年平均气温做减法（-）运算以便获得全球气温变化趋势。

本案例中用到了空间参考投影变换, 栅格分辨率修改, 像素值重分类和栅格代数运算四个功能。

1.2.1 空间参考重投影

ST_Transform函数对栅格数据做空间参考的变换。数据来源不同导致数据的空间参考不同，通过本函数可以将不同来源的数据统一到同一个空间参考系统中。

1.2.2 空间分辨率修改

ST_Resize函数根据用户指定的尺寸和重采样方法对栅格数据进行变换，变化结果对应的地理空间范围保持不变。数据来源的不同导致栅格数据的空间分辨率不尽相同，通过本函数可以将不同的栅格数据确保具备有相同的空间分辨率以便进行下一步计算分析。

1.2.3 像素值重分类

ST_Reclassify函数按照用户指定的规则对栅格数据的像素值进行分类，从而获得一个新的栅格数据。

1.2.4 地图代数运算

ST_MapAlgebra函数使用特定的代数计算表达式对栅格数据的每个像素值进行计算，获得一个新的栅格数据。借助于强大的代数计算表达式，用户可以非常方便地对栅格数据进行运算操作

????2. 实战步骤

以海口市台风“海鸥”(201418)洪水灾害为例，计算承灾体洪水灾害损失率和损失价值分布。具体计算步骤如下：

????2.1 数据入库

数据入库时可以根据需要对数据进行预处理，如使用ST_Transform进行投影变换、ST_Resize修改分辨率等操作，从而获得指定的空间参考和分辨率：
*左右滑动阅览

-- 数据表
create table loss(name varchar(20), rast raster);

-- 导入为指定空间参考和空间分辨率
-- 土地利用类型
INSERT INTO loss values('land_type', ST_Resize(
    ST_Transform(
        ST_CreateRast('<OSS_PATH>/landuse_type_haikou_2015_Albers_30m.tif'), 
            4326, 
      '{"resample":"Near","nodata":true}', 
      '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'), 
    1000, 
    1000, 
    '{"resample":"Near","nodata":true}', 
    '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'));

-- 洪水深度
INSERT INTO loss values('flood_height', ST_Resize(
  ST_Transform(
     ST_CreateRast('<OSS_PATH>/TC_ID_Uniq_201418_Kalmaegi_flood_depth_max_WGS84.tif'), 
        4326, 
        '{"resample":"Near","nodata":true}', 
        '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'), 
     1000, 
     1000, 
     '{"resample":"Near","nodata":true}', 
     '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'));


-- 土地价值
INSERT INTO loss values('land_value', ST_Resize(
  ST_Transform(
     ST_CreateRast('<OSS_PATH>/landuse_value_haikou_2015_Albers_30m.tif'), 
        4326, 
        '{"resample":"Near","nodata":true}', 
        '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'), 
     1000, 
     1000, 
     '{"resample":"Near","nodata":true}', 
     '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'));

????2.2 承灾体损失率计算
1，对不同土地利用进行赋值，属于该类别赋值1，否则为0，获得单种土地利用分布图
*左右滑动阅览

 -- 以type=21 有林地为例
 INSERT INTO loss (name, rast)
 SELECT 'land_type21', ST_Reclassify(rast, '[{"band":0, "remap":{"[0,21)":"0","21":"1","(21,254]":"0"}, "nodata":true, "nodatavalue":255}]', '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}')
 FROM loss
 WHERE name='land_type';

2，洪水淹没水深与各类土地利用赋值结果相乘，得出各类土地利用淹没深度

*左右滑动阅览

-- 使用表达式[0,0] * [1,0] 进行相乘操作
WITH foo as
 (
  SELECT rast from loss WHERE name in ('flood_height', 'land_type21' )
 )
 INSERT INTO loss (name, rast)
 SELECT 'flood_height_type21', ST_MapAlgebra(ARRAY(select rast from foo),
 '[{"expr":"([0,0] * [1,0])","nodata": true, "nodataValue":255}]',
 '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table", "celltype":"64BF"}');

3，按照各类土地利用淹没水深-损失率对照表，对淹没深度进行重分类，赋予对应土地利用类型相应的损失率

*左右滑动阅览

-- 按照洪水深度重新赋值
INSERT INTO loss (name, rast)
SELECT 'ratio_type21', ST_Reclassify(rast, '[{"band":0, "remap":{"[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.9,1.2,1.5,1.8,2.1,2.4,2.7,3,254)":"0,1.5,2.5,4,7,10,13,15,20,24,27,30,30,30,30,30"}, "nodata":true, "nodatavalue":255}]', '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}')
FROM loss
WHERE name='flood_height_type21';

4，有林地土地利用类型对应的损失率图效果

????2.3 承灾体损失值计算

依照公式2，将承灾体价值和损失率进行栅格相乘，得到承灾体损失价值分布：
*左右滑动阅览

WITH foo as
 (
  SELECT rast from loss WHERE name in ('land_value', 'ratio_type21' )
 )
 INSERT INTO loss (name, rast)
 SELECT 'loss_type21', ST_MapAlgebra(ARRAY(select rast from foo), 
 '[{"expr":"([0,0] * [1,0])","nodata": true, "nodataValue":255}]',
 '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table", "celltype":"64BF"}');

有林地洪灾损失价值分布图如下：

????2.4 单SQL执行

当然，也可以把之前所有的SQL都组合到一起直接进行计算，同样以有林地为例：
*左右滑动阅览

-- 把以上几个步骤串联到一起
-- 对结果数据的空间参考和分辨率进行预定义
-- 数据写入临时 tmp_chunk_table 表中，后续可以进行删除
WITH loss AS (
    SELECT 'land_type' as name, ST_Resize(
  ST_Transform(
    ST_CreateRast('<OSS_PATH>/landuse_type_haikou_2015_Albers_30m.tif'), 
      4326, -- 最终结果SRID
      '{"resample":"Near","nodata":true}', 
      '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'), 
    1000,  -- 最终结果分辨率X
    1000, -- 最终结果分辨率Y
    '{"resample":"Near","nodata":true}', 
    '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}') AS rast
UNION ALL
    SELECT 'flood_height' as name, ST_Resize(
  ST_Transform(
     ST_CreateRast('<OSS_PATH>/TC_ID_Uniq_201418_Kalmaegi_flood_depth_max_WGS84.tif'), 
        4326, -- 最终结果SRID
        '{"resample":"Near","nodata":true}', 
        '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'), 
     1000, -- 最终结果分辨率X
     1000, -- 最终结果分辨率Y
     '{"resample":"Near","nodata":true}', 
     '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}') AS rast
UNION ALL
    SELECT 'land_value' as name, ST_Resize(
  ST_Transform(
     ST_CreateRast('<OSS_PATH>/landuse_value_haikou_2015_Albers_30m.tif'), 
        4326, -- 最终结果SRID
        '{"resample":"Near","nodata":true}', 
        '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}'), 
     1000, -- 最终结果分辨率X
     1000, -- 最终结果分辨率Y
     '{"resample":"Near","nodata":true}', 
     '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}') AS rast
),
land_type AS (
SELECT ST_Reclassify(rast, 
              '[{"band":0, "remap":{"[0,21)":"0","21":"1","(21,254]":"0"},"nodata":true, "nodatavalue":255}]', 
          '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}') AS rast
  FROM loss
  WHERE name='land_type'), 
flood_height AS (
SELECT ST_MapAlgebra(ARRAY(select rast from land_type 
                           UNION ALL 
                           select rast from loss WHERE name='flood_height'), 
                 '[{"expr":"([0,0] * [1,0])","nodata": true, "nodataValue":255}]',
                 '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table", "celltype":"64BF"}') AS rast
),
loss_ratio AS (
SELECT ST_Reclassify(rast, 
                        '[{"band":0, "remap":{"[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.9,1.2,1.5,1.8,2.1,2.4,2.7,3,254)":"0,1.5,2.5,4,7,10,13,15,20,24,27,30,30,30,30,30"}, "nodata":true, "nodatavalue":255}]', 
            '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table"}') AS rast
    FROM flood_height),
foo AS (
  SELECT rast from loss WHERE name in ('land_value')
UNION ALL
  SELECT rast from loss_ratio
)
SELECT 'loss_type21', ST_MapAlgebra(ARRAY(select rast from foo), 
    '[{"expr":"([0,0] * [1,0])","nodata": true, "nodataValue":255}]',
    '{"chunking":true,"chunkdim":"(256,256,1)","compression":"lz4","interleaving":"bsq","chunktable":"tmp_chunk_table", "celltype":"64BF"}');

是不是感觉很神奇呢？

????3. 总结

利用Ganos的时空栅格存储、计算和分析能力，将复杂的洪涝承灾体损失计算模型转化为简单的Geo-SQL语句，使得过去必须借助于GIS软件的专业的时空数据处理流程能在数据库内实现，简化用户的程序逻辑，降低开发复杂度与维护成本，使云GIS能力赋能行业用户。

码农公寓