實戰課堂 | 科技減災如何使用Geo-SQL實作洪災承災體損失評估？

🔸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-PG産品。

其中針對時空栅格資料類型，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能力賦能行業使用者。

實戰課堂 | 科技減災如何使用Geo-SQL實作洪災承災體損失評估？

🔸1. 簡介

🔺1.1 模型簡介

1.1.1 資料輸入

1.1.2 資料計算

1.1.3 資料輸出

🔺1.2 Ganos

1.2.1 空間參考重投影

1.2.2 空間分辨率修改

1.2.3 像素值重分類

1.2.4 地圖代數運算

🔸2. 實戰步驟

🔺2.1 資料入庫

🔺2.3 承災體損失值計算

🔺2.4 單SQL執行

🔸3. 總結

繼續閱讀

報錯：'mysql' 不是内部或外部指令，也不是可運作的程式或批處理檔案。

Linxu常用指令技巧彙總

ERROR 1 (HY000): Can't create/write to file '/tmp/#sql_4188_1.MYI' (Errcode: 28)

艱難安裝LDAP,SSL認證

《Linux指令行與Shell腳本程式設計大全第2版.布盧姆》pdf

MySQL的4種隔離級别？出現問題

XX系統實施過程問題總結

無元件上傳圖檔到資料庫中，最完整解決方案

【MySQL資料庫】資料庫索引事務1.索引2.事務

neo4j之cypher使用文檔

【python】【資料處理】畫多元資料分布圖

NOSQL安全攻擊

mybatis_入門程式Mybatis入門

登入plsql 報錯 the account is locked --使用者被鎖

sqlServer根據經緯查距離

SequoiaDB巨杉資料庫C++驅動概述

實戰課堂 | 科技減災 如何使用Geo-SQL實作洪災承災體損失評估？

🔸1. 簡介

🔺1.1 模型簡介

1.1.1 資料輸入

1.1.2 資料計算

1.1.3 資料輸出

🔺1.2 Ganos

1.2.1 空間參考重投影

1.2.2 空間分辨率修改

1.2.3 像素值重分類

1.2.4 地圖代數運算

🔸2. 實戰步驟

🔺2.1 資料入庫

🔺2.3 承災體損失值計算

🔺2.4 單SQL執行

🔸3. 總結

繼續閱讀

實戰課堂 | 科技減災如何使用Geo-SQL實作洪災承災體損失評估？