本文介紹GanosBase的柵格引擎（GanosBase Raster）在水利/氣象領域的分析情境應用。GanosBase通過在資料庫中原生內建影像與格網資料的儲存、檢索與分析能力，為氣象、水利、資源管理、應急及傳媒等領域的提供海量柵格資料的分析與挖掘能力。
關於GanosBase Raster
柵格資料
柵格資料（Raster Data）是一種將現實空間劃分為規則網格的資料形式，每個網格稱為單元（像元或像素），並在各單元上賦予相應的屬性值，以表示實體的特徵。
柵格資料通常有兩種類型：
專題資料：每個柵格像元的值可以是測量值或分類值，例如高程資料（DEM）、汙染物濃度、訊號強度、降雨量、土地權屬類型或植被類型等。
影像資料：又稱遙感影像或遙感映像（Remote Sensing Image），指通過地面遙感、航空遙感或航天遙感平台擷取的映像，這些映像記錄各種地物所反射的電磁傳輸速率征，包括航空影像和遙感衛星影像等。
每一幅柵格資料均具備時間屬性和空間屬性，稱之為時空柵格。時空柵格進一步強調了柵格資料的時態特性，例如對時間序列的系列柵格資料進行管理。
什麼是GanosBase Raster
GanosBase Raster是對象關係型資料庫PostgreSQL相容版本（PolarDB PostgreSQL版（相容Oracle））的一個時空引擎擴充，使上述資料庫能夠有效快速儲存管理柵格類型資料，同時支援多源柵格資料（如遙感、攝影測量和專題地圖）之間的融合與分析，並提供了GeoServer外掛程式協助使用者將庫內柵格對象發布為OGC標準的服務（如WMS或WMTS等）。GanosBase Raster可用於包括氣象、環境監測、地質勘探、自然資源管理、國防、應急響應、電信、傳媒、交通、城市規劃以及國土安全等領域。
模型概述
GanosBase Raster資料模型主要由以下幾個元素構成：
Raster：泛指一份柵格資料。例如，一個景遙感影像、一個TIFF檔案等。
Tile：資料分塊，為一系列像素的集合。Tile為Raster對象在資料庫中儲存的基本單元，一般每個Tile包含256x256個像素值。
Band：由多個Tile構成的2D柵格資料圖層，每個Tile擁有一個行列號。
Cell：表示Tile中的一個像素，可以擁有不同的資料類型。例如Byte，Short，Int，Double等。
Pyramid：通過逐級抽稀的柵格金字塔，方便快速顯示。每個Pyramid包含不同的層級，每個層級對應一個Layer，第0層代表未經處理資料。
Metadata：柵格的儲存中繼資料，包括：空間範圍、投影類型、像素類型等。
如上圖所示，GanosBase Raster採用了一種簡單而高效的通用柵格資料模型來管理專題資料和遙感影像資料。一幅柵格資料（Image）在資料庫中以柵格對象（Raster）形式進行儲存。Raster對象邏輯上由若干可以表示為2D柵格圖層的波段（Band）組成，各個Band的資訊會在入庫過程中從原始影像資料讀取。Raster是以資料區塊（Tile）為基本存放裝置單元進行儲存和管理的，Tile的尺寸預設為256x256像素，但也可以由使用者進行定義。每個Tile可以包含一個或者多個Band。Tile中的一個像素由一個像素單元（Cell）表示。每一個Raster對象都有對應的中繼資料（Metadata）資訊，如圖幅範圍（Extent）、資料類型、投影資訊、行列號等。 如果對柵格資料建立了金字塔模型，則一個Band會包含多個層級（level）的金字塔資料。
功能優勢
相較於開源空間資料庫PostGIS內建的Raster外掛程式，GanosBase Raster在業務適配、儲存成本、計算能力方面都有明顯的優勢，主要體現在以下幾方面：
GanosBase Raster儲存結構更加貼近業務需求。
與PostGIS Raster完全柵格化的儲存方式不同，GanosBase Raster提供了物件導向的儲存結構，一幅柵格資料（影像，映像、DEM等）匯入後直接對應GanosBase Raster的一行記錄，對使用者而言形成了一對一的關係，清晰明了，單個對象的容量沒有限制，一行記錄可以儲存大到1 TB以上的超大影像/柵格。GanosBase Raster屏蔽了直接對Tile進行操作的方式，可完整記錄柵格對象的中繼資料資訊，並且與時序高度關聯，可以與業務進行更好的銜接。
GanosBase Raster支援資料降本儲存。
由於獨特的結構設計，針對海量影像分析這類儲存成本較大的情境，GanosBase Raster支援將柵格中繼資料存放區在庫內，柵格屬性資料存放區在更為廉價的Object Storage Service上。這種情況下依然支援對柵格進行各類空間分析操作，同時可大幅度降低儲存成本。
GanosBase Raster擁有更多時空運算元。
除了支援傳統的柵格空間關係判斷、柵格金字塔、柵格像素值、柵格屬性、柵格影像處理等操作外，GanosBase Raster還支援多種獨特的柵格統計、柵格代數操作，同時支援專業的映像勻色演算法和海量柵格概視圖實現渲染加速。
面雨量分析應用情境
情境描述
面雨量是描述整個地區（流域）內單位面積上平均降水量的物理量，能夠較為客觀地反映該地區的降水情況。面雨量是洪水預報中極為重要的指標，為各級政府在防汛抗洪及水庫調蓄方面提供了依據，同時也為防災抗災和經濟建設提供參考。由於面雨量作為地區性指標，相關管理部門通常只能擷取觀測網站提供的點位監測資料，因此需要構建從點位元據到面指標的計算鏈路，並提供相關附加分析資料以輔助決策。
需求分析
為實現上述情境，需要考慮基於離散監測點進行空間插值，以構建格網資料。隨後，基於該網格資料開展等值線/面追蹤，同時計算任意向量面所覆蓋範圍內的相關統計資訊，並求取總體面雨量。
資料匯入：使用FDW擴充實現向量資料的入庫，完成向量觀測點資料的入庫。
資料插值：使用ST_InterpolateRaster函數基於向量觀測點插值為自訂大小格網，用於儲存GanosBase Raster。
等值線/面：使用ST_Contour函數基於上述格網資料追蹤等值線或面，以分析總體雨量的分布情況。
空間統計：根據不同格網中的雨量分布進行分類，使用ST_Statistics函數統計各類別的最大雨量、最小雨量、中位元及總雨量等資訊。
面雨量計算：基於任意向量邊界與格網，計算出向量所包含的格網。針對穿透的格網，構建相應的統計規則，最終得出該向量面的總體雨量資料。
技術實現
安裝GanosBase外掛程式
在目標資料庫中安裝ganos_raster和ganos_fdw外掛程式。
CREATE EXTENSION ganos_raster CASCADE;
CREATE EXTENSION ganos_fdw CASCADE;匯入資料
準備shapefile格式的點要素圖層檔案point.shp和面要素圖層檔案polygon.shp，其中每個點要素具有一個id（字串）屬性和一個降雨量pre（浮點型）屬性。將圖層疊加至QGIS中顯示，如下圖所示，數字表示每個點的監雨量屬性值：
通過GanosBase FDW將shapefile檔案以外表的形式從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>'); -- 將shapefile映射為外表 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;通過外表映射的方式匯入資料後，可以使用SQL進行各類查詢，例如查詢point_table中的所有記錄：
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，支援將點要素插值為Raster類型的對象。空間插值函數同樣支援並行操作，能夠通過配置並行度來提升執行效率。
建立包含有Raster欄位類型的表，用於儲存插值後的柵格資料。
CREATE TABLE IF NOT EXISTS raster_table ( id integer, rast raster -- raster資料類型，用於儲存插值後的柵格對象 );調用ST_InterpolateRaster函數對point_table表中的空間點資料採用反距離加權（IDW）方法進行插值，並將結果插入至raster_table中。
說明樣本中使用ST_MakePoint方法建立一個包含屬性值（降雨量）的Point對象，以作為插值的輸入參數。
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函數將產生的raster對象匯出為檔案形式以便查看，此處採用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中的行政區劃Polygon對象對插值後的柵格對象進行裁剪，產生新的柵格對象，並將其插入到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產生等值線，並與行政區劃Polygon進行求交，最終結果儲存在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屬性並將其指定為true，則ST_Contour將以等值面形式返回結果，具體的SQL語句如下：
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函數將輸出POLYGON類型的面要素，經過在QGIS中疊加配色後，效果如下：
雨量分析統計
以下將展示如何利用前述通過降雨量觀測值插值獲得的柵格對象，計算手繪多邊形的面雨量。面雨量是指在特定地區內，降雨總量與該地區面積的比值，通常以毫米/平方米為單位表示。因此，首先計算每個像素覆蓋的多邊形地區的面積，隨後將該像素值乘以相應面積並進行求和，以獲得總降雨量。最後，將總降雨量除以多邊形的面積即可得到面雨量。
統計任意地區降雨量統計資訊。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.84mm，最小降雨量為2.5mm，而降雨量在大於10mm且小於等於15mm的像素總數為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;返回結果如下：
row | col | st_astext | pre -----+-----+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+-------------------- 107 | 33 | POLYGON((119.244752593338 28.0907410103828,119.244752593338 28.0878755431622,119.2479422912 28.0878755431622,119.2479422912 28.0907410103828,119.244752593338 28.0907410103828)) | 5.66606330871582 107 | 34 | POLYGON((119.244752593338 28.0878755431622,119.244752593338 28.0850100759417,119.2479422912 28.0850100759417,119.2479422912 28.0878755431622,119.244752593338 28.0878755431622)) | 5.698885917663574 107 | 35 | POLYGON((119.244752593338 28.0850100759417,119.244752593338 28.0821446087211,119.2479422912 28.0821446087211,119.2479422912 28.0850100759417,119.244752593338 28.0850100759417)) | 5.7316670417785645 107 | 36 | POLYGON((119.244752593338 28.0821446087211,119.244752593338 28.0792791415006,119.2479422912 28.0792791415006,119.2479422912 28.0821446087211,119.244752593338 28.0821446087211)) | 5.764394283294678 107 | 37 | POLYGON((119.244752593338 28.0792791415006,119.244752593338 28.0764136742801,119.2479422912 28.0764136742801,119.2479422912 28.0792791415006,119.244752593338 28.0792791415006)) | 5.797055244445801 107 | 38 | POLYGON((119.244752593338 28.0764136742801,119.244752593338 28.0735482070595,119.2479422912 28.0735482070595,119.2479422912 28.0764136742801,119.244752593338 28.0764136742801)) | 5.829642295837402 107 | 39 | POLYGON((119.244752593338 28.0735482070595,119.244752593338 28.070682739839,119.2479422912 28.070682739839,119.2479422912 28.0735482070595,119.244752593338 28.0735482070595)) | 5.862148761749268 107 | 40 | POLYGON((119.244752593338 28.070682739839,119.244752593338 28.0678172726184,119.2479422912 28.0678172726184,119.2479422912 28.070682739839,119.244752593338 28.070682739839)) | 5.894571304321289 107 | 41 | POLYGON((119.244752593338 28.0678172726184,119.244752593338 28.0649518053979,119.2479422912 28.0649518053979,119.2479422912 28.0678172726184,119.244752593338 28.0678172726184)) | 5.92690896987915 107 | 42 | POLYGON((119.244752593338 28.0649518053979,119.244752593338 28.0620863381773,119.2479422912 28.0620863381773,119.2479422912 28.0649518053979,119.244752593338 28.0649518053979)) | 5.959161758422852 (10 rows)計算面雨量。面雨量單位按照每平方米進行計算，因此在計算每個像素多邊形的面積之前，需要通過ST_Transform方法將座標系從EPSG:4326轉換為EPSG:3857，以獲得以平方米為單位的面積值。然後，將該像素的像素值相乘，最後對所有像素的結果進行求和，併除以所有像素的總面積，以計算出面雨量。具體SQL語句如下：
SELECT sum(ST_Area(ST_Transform(geom,3857)) * pre) / sum(ST_Area(ST_Transform(geom,3857))) as "precipitation(mm/m^2)" FROM pixel_pre WHERE pre>0;輸出面雨量結果如下：
precipitation(mm/m^2) ----------------------- 9.434021577002174 (1 row)
總結
GanosBase柵格引擎提供了一整套針對格網資料的匯入、儲存、分析與可視化支撐能力，協助高效解決遙感影像管理、DEM資料分析、格網資料分析等業務，同時還具備GPU運算等高階能力，以支援不同的應用需求。GanosBase柵格引擎在面雨量計算情境中的應用，已有效支撐水利部氣象降水會商系統的運行，並將該系統的分析計算效率提升了十倍以上。目前，GanosBase柵格引擎在水利、自然資源、氣象、環境保護、應急等領域具備全面的能力供給與應用案例，能夠為多家客戶的空天巨量資料管理類應用提供穩定、高效且健全的時空基底保障。