GanosBase ラスターエンジンに基づく面積雨量解析 - PolarDB

このトピックでは、水利気象解析において GanosBase（GanosBase Raster とも呼ばれます）のラスターエンジンを使用する方法について説明します。GanosBase は、データベース内での画像データとグリッドデータの格納、取得、および解析を統合し、気象、水利、資源管理、緊急対応、メディアなどの分野で大規模なラスターデータの解析と処理を可能にします。

概要GanosBase Raster について

ラスターデータ

ラスターデータは、現実世界の空間をグリッドに分割します。各グリッドはセルまたはピクセルと呼ばれ、エンティティの特性を表すプロパティが含まれています。

ラスターデータは、通常、次の 2 つのカテゴリに分類されます。

主題データ: 各ラスターピクセルの値は、数値標高モデル（DEM）データ、汚染物質濃度、電波強度、降水量、土地所有権、植生タイプなどの情報を記述するための測定値または分類です。
画像データ（リモートセンシング画像とも呼ばれます）: 画像は、地上、航空、または衛星のリモートセンシングプラットフォームによってキャプチャされます。これらの画像は、さまざまなオブジェクトの電磁的特性を反映しており、航空写真と衛星写真の両方を含みます。

ラスターデータには空間属性と時間属性の両方が含まれているため、ラスターは時空間ラスターとも呼ばれます。時間属性を使用すると、ラスターデータを使用して時系列を管理できます。

GanosBase Raster とは

GanosBase Raster は、PolarDB for PostgreSQL によって提供される時空間エンジン拡張機能です。この拡張機能により、データベースはラスターデータ型を効率的に格納、管理、および処理できます。リモートセンシング、写真測量、主題マッピングなど、さまざまなソースからのラスターデータの統合と解析をサポートしています。また、ラスターオブジェクトを WMS や WMTS などの標準の Open Geospatial Consortium（OGC）サービスとして公開するための GeoServer 拡張機能も提供しています。GanosBase Raster は、気象、環境モニタリング、地質探査、天然資源管理、国防、緊急対応、電気通信、メディア、輸送、都市計画、国土安全保障など、幅広い分野で使用できます。

モデルの概要

GanosBase ラスターデータモデルは、次の要素で構成されています。

ラスター: リモートセンシング画像や TIFF ファイルなどのラスターデータセット。
タイル: ピクセルのブロック。タイルは、ラスターオブジェクトの基本的なストレージ単位です。各タイルは、デフォルトで 256 x 256 ピクセルです。
バンド: ラスターデータセット内の単一の 2D レイヤー。各バンドは複数のタイルで構成され、各タイルには独自の座標番号があります。
セル: タイル内のピクセル。セルは、バイト、short、int、double など、さまざまなデータ型にすることができます。
ピラミッド: ラスターデータの表示を高速化するために使用される、ラスターの解像度を下げた一連のバージョン。ピラミッドにはさまざまなレベルがあり、各レベルはレイヤーに対応しています。レベル 0 は生データを参照します。
メタデータ: 空間範囲、投影タイプ、ピクセルタイプなどのラスターメタデータ。

前の図に示すように、GanosBase Raster は、シンプルながらも効率的なラスターモデルを使用して、主題データとリモートセンシング画像データを管理します。データベースシステムでは、ラスター画像はラスターオブジェクトとして格納されます。ラスターオブジェクトは論理的に複数のバンドに分割され、各バンドはインポートプロセス中に取得された生の画像データを保持します。ラスターオブジェクトは、タイルによって格納および管理されます。タイルのサイズはデフォルトで 256 x 256 ピクセルですが、必要に応じて変更できます。各タイルには、1 つ以上のバンドが含まれています。タイルは、ピクセルを表すセルで構成されています。各ラスターオブジェクトには、空間範囲、データ型、投影情報、座標番号などの独自のメタデータがあります。ラスターがピラミッド構造で編成されている場合、各バンドには独自のピラミッドがあります。

メリット

GanosBase Raster は、ビジネス適応、ストレージの節約、計算機能の点で、PostGIS Raster よりも優れています。

優れたビジネス適応性
PostGIS Raster はデータを完全にラスター化しますが、GanosBase Raster はオブジェクト指向のストレージ構造を提供します。画像でも DEM ファイルでも、各ラスターは単一のデータベースレコードとして格納され、1 対 1 のマッピングを形成します。これらの各レコードには、最大 1 TB のラスターを含めることができます。さらに、GanosBase ではタイルに対する直接操作は許可されていないため、メタデータの整合性が保証されます。さらに、時系列データと高度に相関しています。つまり、GanosBase Raster を使用すると、ラスターデータをビジネスプロセスにより適切に統合できます。
ストレージの大幅な節約
GanosBase Raster は、ラスターメタデータをデータベースに、ラスター属性データを低コストの OSS（Object Storage Service）に保存するアーキテクチャを採用しています。これは、大量の画像の解析など、大容量のストレージを必要とするユースケースに役立ちます。このようにして、ストレージコストを大幅に削減しながら、ラスターの空間解析を実行できます。
より多くの時空間演算子
空間関係の評価、画像の処理、ピラミッド、ピクセル値、属性に対する標準操作に加えて、GanosBase Raster は、さまざまな独自の統計演算と代数演算、特殊な画像カラーバランスアルゴリズム、大規模なラスターデータセットを処理する場合でも概要マップの高速レンダリングを提供します。

面積雨量解析

シナリオの説明

面積雨量は、洪水予測の重要な指標であり、地域内の単位面積あたりの平均降水量を反映しています。洪水対策、貯水池の調整、防災、経済発展のための重要なデータを提供します。ほとんどの場合、関連部門は孤立した観測地点からモニタリングデータを収集するため、意思決定のために地域全体の面積雨量解析を実行するための計算リンクを確立する必要があります。

要件

このシナリオでは、グリッドデータを作成するために、離散モニタリングポイントの空間補間が必要です。グリッドデータに基づいて、等高線とサーフェストレーシングが実行されます。同時に、ベクトルサーフェスで覆われた範囲内の関連する統計情報を計算して、全体的な面積雨量を決定します。

データのインポート: FDW 拡張機能を使用して、観測ポイントのベクトルデータをデータベースにインポートします。
データの補間: ST_InterpolateRaster 関数を使用して、ベクトル観測ポイントを指定されたサイズのグリッドに変換します。これは、GanosBase Raster データを格納するために使用されます。
等高線とサーフェス: ST_Contour 関数を使用して、グリッドデータに基づいて等高線またはサーフェスをトレースし、全体的な降水分布を解析します。
空間統計: ST_Statistics 関数を使用して、グリッド全体の降水分布を分類し、各カテゴリの最大、最小、中央値、および合計降水量を決定します。
面積雨量の計算: 任意のベクトル境界内のグリッドを計算し、貫通グリッドの統計ルールを作成して、ベクトル領域の全体的な降水量データを取得します。

実装

ack-onepilot コンポーネントをインストールするGanosBase 拡張機能をインストールする拡張機能

ターゲットデータベースに ganos_raster 拡張機能と ganos_fdw 拡張機能をインストールします。

CREATE EXTENSION ganos_raster CASCADE;
CREATE EXTENSION ganos_fdw CASCADE;

データのインポート

point.shp という名前のポイントレイヤーファイルと polygon.shp という名前のポリゴンレイヤーファイルを準備します。どちらもシェープファイル形式です。各ポイントには、ID 属性（文字列）と降水量属性（浮動小数点数）が含まれています。次の例に示すように、QGIS に重ねて表示されたレイヤーを表示します。各数値は、各ポイントの降水量を示しています。

GanosBase FDW を使用して、シェープファイルを OSS からデータベースへの外部テーブルとしてマップします。次に、CREATE 文を実行してテーブルを作成し、データを挿入します。詳細については、「FDW」をご参照ください。次の SQL ステートメントの例を参照してください。

-- ak と ak_secret は、OSS サービスの AccessKeyId と AccessKeySecret です。
CREATE SERVER myserver
  FOREIGN DATA WRAPPER ganos_fdw
  OPTIONS (
    datasource 'OSS://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file',
    format 'ESRI Shapefile' );
CREATE USER MAPPING FOR CURRENT_USER SERVER myserver OPTIONS (user '<ak>', password '<ak_secret>');

-- シェープファイルを外部テーブルとしてマップする
CREATE FOREIGN TABLE foreign_point_table (
  fid integer,
  id  varchar,
  geom geometry,
  pre double precision -- 降水量)
  SERVER myserver
  OPTIONS ( layer 'point' );

CREATE FOREIGN TABLE foreign_polygon_table (
  fid integer,
  id  varchar,
  geom geometry)
  SERVER myserver
  OPTIONS ( layer 'polygon' );

-- テーブルを作成し、データを挿入する
CREATE TABLE point_table AS 
  SELECT fid, geom, pre FROM foreign_point_table;

CREATE TABLE polygon_table AS 
  SELECT fid, geom FROM foreign_polygon_table;

外部テーブルを介してデータをインポートした後、point_table からすべてのレコードをクエリするなど、さまざまな SQL クエリを実行できます。

SELECT fid, ST_AsText(geom),pre FROM point_table;

結果の例:

fid  |          st_astext          | pre  
-----+-----------------------------+------
   0 | POINT(119.1084 28.50302)    |    5
   1 | POINT(118.768925 28.475747) |  3.5
   2 | POINT(119.30954 28.773729)  |  2.5
   3 | POINT(119.039694 28.363413) |    6
   4 | POINT(119.035561 28.614094) |    4
   5 | POINT(119.9517 28.77)       |  0.5
   6 | POINT(120.35833 28.62694)   |    1
   7 | POINT(119.908078 28.439481) |  1.5
   8 | POINT(119.472 28.5933)      |  4.5
   9 | POINT(119.400282 28.398895) |    4
  10 | POINT(119.783954 28.271403) |    0
  11 | POINT(119.663102 28.514025) |    3
  12 | POINT(120.343889 28.9031)   |    0
  13 | POINT(119.425841 27.768324) |  8.5
  14 | POINT(119.426237 28.015453) |  6.5
  15 | POINT(119.677214 27.944789) |  4.5
  16 | POINT(119.078999 28.045044) |    7
  17 | POINT(120.328711 28.411951) |  1.5
  18 | POINT(120.071226 28.412596) |    1
  19 | POINT(120.336292 27.986628) |    2
  20 | POINT(119.174307 27.635898) |   17
  21 | POINT(119.106197 27.776089) | 13.5
  22 | POINT(119.316833 28.191166) |    4
  23 | POINT(119.554612 28.149524) |  3.5
  24 | POINT(119.573099 28.295702) |  2.5
  25 | POINT(119.154161 28.237417) |    5
  26 | POINT(120.391447 28.211369) |    4
  27 | POINT(119.979144 28.568005) |    1
  28 | POINT(119.468737 27.606996) | 10.5
  29 | POINT(119.881658 28.03036)  |  1.5
  30 | POINT(120.068508 28.165111) |  1.5
  31 | POINT(119.765358 27.821914) |    2
  32 | POINT(118.892389 27.733222) |   14
  33 | POINT(118.75833 28.00694)   |  4.5
  34 | POINT(118.91 27.51222)      | 11.5
  35 | POINT(119.248086 27.435982) | 14.5
(36 rows)

空間補間を実行する

GanosBase は、空間補間関数 ST_InterpolateRaster を提供します。これにより、補間によってポイントをラスターオブジェクトに変換できます。この関数は、並列処理の次数を設定することで、実行パフォーマンスを向上させる並列処理をサポートしています。

補間されたラスターデータを格納するためのラスターフィールドを持つテーブルを作成します。

CREATE TABLE IF NOT EXISTS raster_table
(
  id integer, 
  rast raster -- 補間されたラスターオブジェクトを格納するためのラスターデータ型
);

ST_InterpolateRaster 関数を使用して、逆距離加重（IDW）メソッドを使用して point_table から空間データを補間し、結果を raster_table に挿入します。
説明
この例では、ST_MakePoint を使用して、降水量などの属性を含むポイントオブジェクトを生成し、補間の入力パラメーターとして使用します。
```
INSERT INTO raster_table(id, rast) VALUES(
  1, 
  (SELECT ST_InterpolateRaster(
    ST_Collect(ST_MakePoint(ST_X(geom),ST_Y(geom),pre)),
    512,
    512,
    '{"method":"IDW","radius":2.0,"max_points":"30","min_points":"1"}',
    '{"chunktable":"rbt","celltype":"32bf"}')
  FROM point_table));
```
大量のポイントを処理する場合は、補間関数を実行する前に並列処理を有効にして、パフォーマンスを向上させます。
```
-- 並列処理の次数を 4 に設定する
SET ganos.parallel.degree = 4;
```
ST_ExportTo 関数を使用して、生成されたラスターオブジェクトをファイルにエクスポートして表示します。この場合、ファイルは COG 形式で OSS にエクスポートされます。
```
SELECT ST_ExportTo(rast, 
  'COG',
  'OSS://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file.tif') 
  FROM raster_table WHERE id=1;
```
ファイルをダウンロードし、QGIS で開いて結果を表示します。
ST_ClipToRast 関数を使用して、補間されたラスターオブジェクトを polygon_table のポリゴンオブジェクト（行政区画を表す）でクリップします。このプロセスでは、新しいラスターオブジェクトが作成され、オブジェクトが raster_table に挿入されます。
```
INSERT INTO raster_table
SELECT 2,ST_ClipToRast(r.rast, p.geom, 0, '', NULL, '', '{"chunktable":"rbt"}') 
  FROM raster_table AS R, polygon_table AS p;
```
次のステートメントを実行して、クリップされたラスターデータを idw_clip.tif という名前のファイルにエクスポートします。
```
SELECT ST_ExportTo(rast, 
  'COG',
  'OSS://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file.tif') 
  FROM raster_table WHERE id=2;
```

輪郭とサーフェスのトラッキング

補間によってラスターオブジェクトが生成された後、GanosBase が提供するラスター処理関数を使用して、データの解析と処理を行うことができます。次の例は、等高線とサーフェストレーシングを実行する方法を示しています。

等高線を生成する

GanosBase は、ST_Contour 関数を提供して、等高線とサーフェスを作成します。次の SQL 文を使用して、1.0 の等間隔で等高線を生成し、それらを行政区画ポリゴンに重ね合わせます。結果は rs_contours テーブルに格納されます。

CREATE TABLE rs_contours AS
SELECT id, max_value, min_value, ST_Intersection(a.geom,p.geom) FROM
(SELECT (ST_Contour(rast,1,'{"interval":"1.0"}')).*
FROM raster_table WHERE id=1) a, polygon_table AS p;

作成した等高線テーブルである rs_contours を QGIS にインポートして、オーバーレイ表示を行います。次の結果の画像では、緑色の数字は観測サイトのメジャーであり、赤色の数字は等高線の値です。

サーフェスを生成する

polygonize 属性を ST_Contour 関数に追加し、true に設定してサーフェスを返します。

CREATE TABLE rs_contours_polygon AS
SELECT id, max_value, min_value, ST_Intersection(a.geom, p.geom) FROM
(SELECT (ST_Contour(rast, 1, '{"interval": "1.0", "polygonize": "true"}')).*
FROM raster_table WHERE id = 1) a, polygon_table AS p;

ST_Contour 関数は、ポリゴンタイプのサーフェスを出力します。QGIS でオーバーレイされ、色付けされたサーフェスの次のサンプルイメージをご覧ください。

降水量分析

このセクションでは、補間されたラスタオブジェクトを使用して、特定のポリゴンの面積降水量を計算するプロセスについて説明します。面積降水量とは、特定のエリア内の総降水量とそのリージョンの面積の比率を指し、通常は平方メートルあたりのミリメートル（mm/m²）で表されます。したがって、まず各ピクセルで覆われたポリゴンの面積が計算されます。次に、ピクセル値に対応する面積を乗算し、合計して総降水量が得られます。最後に、総降水量をポリゴンの面積で割って、面積降水量が計算されます。

エリアの降水量の統計を収集します。GanosBase は、ラスターデータ統計を収集するための ST_Statistics 関数を提供します。この関数は、任意のエリア（ポイント、ライン、およびエリア）内のラスター画素値をサポートしています。任意のポリゴンエリアの出力と、補間されたラスターデータに対して間隔ベースの統計を実行するための画素値統計範囲を指定できます。この例では、統計範囲は (0, 5, 10, 15, 20] です。

SELECT (ST_Statistics(rast, ST_GeomFromText('POLYGON((119.0969 28.0519, 118.9058 27.8942, 119.0502 27.5649, 119.3347 27.6292, 119.4262 27.8775, 119.4927 28.1823, 119.3812 28.1186, 119.0969 28.0519))'), 
  0, '(0, 5, 10, 15, 20]', false)).*
FROM raster_table
WHERE id = 1;

結果例：

name     | band |        min         |        max         |        mean        |        sum         | count |        std         |       median       |        mode        
---------+------+--------------------+--------------------+--------------------+--------------------+-------+--------------------+--------------------+--------------------
 full    |    0 | 2.5009584426879883 | 16.841625213623047 |  6.618387373747887 | 243748.58858776093 | 36829 | 2.6857099819905033 |  6.045845985412598 |  3.892089605331421
 (0-5]   |    0 | 2.5009584426879883 |  4.999907970428467 |  4.144016098134749 | 50963.109974861145 | 12298 | 0.5534420165252167 |  4.231814861297607 |  3.892089605331421
 (5-10]  |    0 |  5.000039577484131 |  9.999935150146484 |  6.852788502446667 |  136110.0852355957 | 19862 | 1.2932534814382863 |  6.623260498046875 |  5.028009414672852
 (10-15] |    0 |  10.00007438659668 | 14.993864059448242 | 11.947531793007881 |  52999.25103378296 |  4436 | 1.2279701295194279 | 11.900737762451172 | 12.430845260620117
 (15-20] |    0 |  15.00446891784668 | 16.841625213623047 | 15.777434950734413 |  3676.142343521118 |   233 | 0.5092360295014834 |   15.7352294921875 |  15.00446891784668
(5 rows)

出力された統計から、指定されたポリゴンエリア内の詳細な統計データを取得できます。たとえば、このエリアの最大降水量は 16.84 mm、最小降水量は 2.5 mm で、降水量が 10 mm を超え 15 mm 以下の画素は 4436 個あります。

ST_ClipToRast 関数を使用して、ラスターオブジェクトをクリップし、新しいラスターオブジェクトを生成し、そのオブジェクトを raster_table テーブルに挿入します。
```
INSERT INTO raster_table
SELECT 3, ST_ClipToRast(rast, ST_GeomFromText('POLYGON((119.0969 28.0519, 118.9058 27.8942, 119.0502 27.5649, 119.3347 27.6292, 119.4262 27.8775, 119.4927 28.1823, 119.3812 28.1186, 119.0969 28.0519))'),
  0, '', NULL, '', '{"chunktable":"rbt"}') 
FROM raster_table WHERE id = 1;
```
最も正確なエリア降水量統計を取得するには、ラスターデータの各画素を走査して、画素がカバーするエリアの矩形フィーチャとそれに対応する画素値を決定します。ST_Width 関数と ST_Height 関数を使用して、クリップされたラスターオブジェクトのディメンションを取得します。次の SQL 文を参照してください。
```
SELECT ST_Width(rast), ST_Height(rast) FROM raster_table WHERE id = 3;
```
結果例：
```
st_width  | st_height 
----------+-----------
      185 |       217
```

各画素の面積と降水量を取得します。ST_PixelAsPolygon 関数と ST_Value 関数を使用して、各画素がカバーするポリゴンエリアとそれに対応する画素値を取得し、結果を新しいテーブルに格納します。次の SQL 文を参照してください。

-- 新しいテーブルを作成する
CREATE TABLE pixel_pre
(   
   row integer,
   col integer,
   geom geometry, -- 画素エリア
   pre  double precision -- 画素値
);


UPDATE raster_table SET rast = ST_SetSrid(rast, 4326);

-- クリップされたラスターデータの各画素を走査して、カバーされている矩形エリアと対応する画素値を取得します
DO
$do$
BEGIN 
   FOR i IN 0..184 LOOP
      FOR j IN 0..216 LOOP
            INSERT INTO pixel_pre
            SELECT i, j, ST_PixelAsPolygon(rast, j, i), ST_Value(rast, 0, i, j) as value FROM raster_table where id = 3;
      END LOOP;
   END LOOP;
END
$do$;

新しいテーブルをクエリします。各画素の座標、カバーされている矩形エリア、および値が返されます。

SELECT ROW, col, ST_AsText(geom), pre FROM pixel_pre LIMIT 10;