使用Ganos实现洪涝灾害承灾体损失综合评估_云原生数据库 PolarDB(PolarDB)-阿里云帮助中心

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

简介

利用洪涝灾害承灾体损失综合评估模型，对灾害损失率及损失价值的分布进行科学计算，对于指导洪涝救灾、建立灾害预警机制、深入研究洪涝灾害成灾机制，以及构建和完善更为科学和准确的洪灾损失评估预测体系，具有重要意义。

模型简介

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

数据输入

输入数据包括：

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

数据计算

承灾体损失率
计算公式：DR=f(H) 其中，DR为每类承灾体的损失率，H为致灾因子强度，f为致灾因子和损失率映射关系。
承灾体损失价值
计算公式：Loss=DR*E
其中，Loss为损失价值，DR为灾害损失率，E为承灾体价值。

数据输出

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

Ganos Raster

针对时空栅格数据类型，Ganos提供了超大规格的栅格数据存储与计算能力。单幅栅格数据在理论上不存在容量限制，具备全球范围内的统一管理能力。这使得传统GIS中复杂的栅格分析操作能够通过Geo-SQL轻松实现，同时具备与几何数据类型的一致分析能力，详细介绍可参见栅格模型。

地图代数是栅格分析与地理信息系统（GIS）建模中常用的技术方法。Ganos为栅格图层计算操作提供了栅格代数表达式语言及一系列栅格代数函数，统称为ACL（Algebra Computing Language）。ACL包含通用算术运算符、逻辑运算符、位运算、关系运算符以及一系列统计函数，并允许这些运算符自由组合，以实现更为复杂的运算操作。Ganos栅格利用ACL强大的计算表达式，支持基于一个或多个栅格对象的像元值进行条件查询、数学建模、分类操作以及生成新的结果栅格对象。PL/pgSQL与ACL的结合使用提供更为强大且易用的栅格分析工具。PL/pgSQL支持变量与常量的声明、通用数学表达式、基本函数、逻辑判断及流程控制，而ACL则为栅格计算提供了像元代数计算的表达方式。您可以轻松结合两种优势进行时空栅格的分析与建模，例如对全球年平均气温进行减法运算，以获取全球气温变化趋势。本案例中使用了空间参考投影变换、栅格分辨率修改、像素值重分类和栅格代数运算四个功能。

空间参考重投影

ST_Transform函数用于对栅格数据进行空间参考的变换。由于数据来源的不同，导致其空间参考也有所差异。通过本函数，可以将不同来源的数据统一转换到同一空间参考系统中。详情请参考ST_Transform。

空间分辨率修改

ST_Resize函数依据自定义尺寸和重采样方法对栅格数据进行变换，而变换结果所对应的地理空间范围保持不变。由于数据来源的差异，栅格数据的空间分辨率可能存在不一致性。通过本函数，可以确保不同的栅格数据具备相同的空间分辨率，从而为后续的计算分析奠定基础。详情请参考ST_Resize。

像素值重分类

ST_Reclassify函数根据自定义规则对栅格数据的像素值进行重新分类，以生成一个新的栅格数据。详情请参考ST_Reclassify。

地图代数运算

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

最佳实践

以台风“海鸥”（201418）引发的洪水灾害为例，计算承灾体的洪水灾害损失率及损失价值的分布情况。

数据入库

在数据入库过程中，可以根据具体需求对数据进行预处理，例如采用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://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file1.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://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file2.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://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file3.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"}'));

承灾体损失率计算

对不同土地利用进行赋值：若属于该类别，则赋值为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';

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

-- 使用表达式[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"}');

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

-- 按照洪水深度重新赋值
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';

林地土地利用类型对应的损失率图效果图展示。

承灾体损失值计算

使用承灾体损失价值公式，将承灾体的价值与损失率进行栅格相乘，以获得承灾体损失价值的分布。

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"}');

林地土地利用类型的洪灾损失价值分布图如下：

单SQL执行（可选）

当然，也可以将之前所有的SQL语句组合在一起进行计算，以下以有林地为例：

-- 把以上几个步骤串联到一起
-- 对结果数据的空间参考和分辨率进行预定义
-- 数据写入临时 tmp_chunk_table 表中，后续可以进行删除
WITH loss AS (
    SELECT 'land_type' as name, ST_Resize(
  ST_Transform(
    ST_CreateRast('OSS://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file1.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://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file2.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://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file3.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"}');

总结

Ganos的时空栅格存储、计算与分析能力，能够将复杂的洪涝承灾体损失计算模型转化为简单的Geo-SQL语句。这一创新使得以往需要依赖GIS软件进行的专业时空数据处理流程，能够在数据库内部实现，从而简化程序逻辑，降低开发复杂性与维护成本，进而使云GIS能力更有效地赋能行业用户。