使用GanosBase實現洪澇災害承災體損失綜合評估 - PolarDB

利用GanosBase的時空柵格儲存、計算和分析能力，可將複雜的洪澇承災體損失計算模型轉化為簡單的Geo-SQL語句，使得過去必須藉助GIS軟體的專業的時空資料處理流程能在資料庫內實現，簡化程式邏輯，降低開發複雜度與維護成本，使雲GIS能力賦能行業使用者。

簡介

利用洪澇災害承災體損失綜合評估模型，對災害損失率及損失價值的分布進行科學計算，對於指導洪澇救災、建立災害預警機制、深入研究洪澇災害成災機制，以及構建和完善更為科學和準確的洪災損失評估預測體系，具有重要意義。

模型簡介

洪澇災害承災體損失綜合評估模型利用洪水災害淹沒水深分布，結合承災體類型、承災體價值及脆弱性資料，計算災害損失率和損失價值分布。

資料輸入

輸入資料包括：

洪水淹沒水深分布。
土地利用類型分布。
土地利用價值分布。
土地利用淹沒水深-損失率對照表。

資料計算

承災體損失率
計算公式：DR=f(H) 其中，DR為每類承災體的損失率，H為致災因子強度，f為致災因子和損失率映射關係。
承災體損失價值
計算公式：Loss=DR*E
其中，Loss為損失價值，DR為災害損失率，E為承災體價值。

資料輸出

洪澇災害損失率分布。
洪澇災害損失價值分布。

GanosBase Raster

針對時空柵格資料類型，GanosBase提供了超大規格的柵格資料存放區與計算能力。單幅柵格資料在理論上不存在容量限制，具備全球範圍內的統一管理能力。這使得傳統GIS中複雜的柵格分析操作能夠通過Geo-SQL輕鬆實現，同時具備與幾何資料類型的一致分析能力，詳細介紹可參見柵格模型。

地圖代數是柵格分析與地理資訊系統（GIS）建模中常用的技術方法。GanosBase為柵格圖層計算操作提供了柵格代數運算式語言及一系列柵格代數函數，統稱為ACL（Algebra Computing Language）。ACL包含通用算術運算子、邏輯運算子、位元運算、關係運算子以及一系列統計函數，並允許這些運算子自由組合，以實現更為複雜的運算操作。GanosBase柵格利用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"}');

總結

GanosBase的時空柵格儲存、計算與分析能力，能夠將複雜的洪澇承災體損失計算模型轉化為簡單的Geo-SQL語句。這一創新使得以往需要依賴GIS軟體進行的專業時空資料處理流程，能夠在資料庫內部實現，從而簡化程式邏輯，降低開發複雜性與維護成本，進而使雲GIS能力更有效地賦能行業使用者。