本文介绍Ganos的栅格引擎（Ganos Raster）在水利/气象领域的分析场景应用。Ganos通过在数据库中原生内置影像与格网数据的存储、检索与分析能力，为气象、水利、资源管理、应急及传媒等领域的提供海量栅格数据的分析与挖掘能力。
关于Ganos Raster
栅格数据
栅格数据（Raster Data）是一种将现实空间划分为规则网格的数据形式，每个网格称为单元（像元或像素），并在各单元上赋予相应的属性值，以表示实体的特征。
栅格数据通常有两种类型：
专题数据：每个栅格像元的值可以是测量值或分类值，例如高程数据（DEM）、污染物浓度、信号强度、降雨量、土地权属类型或植被类型等。
影像数据：又称遥感影像或遥感图像（Remote Sensing Image），指通过地面遥感、航空遥感或航天遥感平台获取的图像，这些图像记录各种地物所反射的电磁波特征，包括航空影像和遥感卫星影像等。
每一幅栅格数据均具备时间属性和空间属性，称之为时空栅格。时空栅格进一步强调了栅格数据的时态特性，例如对时间序列的系列栅格数据进行管理。
什么是Ganos Raster
Ganos Raster是对象关系型数据库PostgreSQL兼容版本（PolarDB PostgreSQL版（兼容Oracle））的一个时空引擎扩展，使上述数据库能够有效快速存储管理栅格类型数据，同时支持多源栅格数据（如遥感、摄影测量和专题地图）之间的融合与分析，并提供了GeoServer插件帮助用户将库内栅格对象发布为OGC标准的服务（如WMS或WMTS等）。Ganos Raster可用于包括气象、环境监测、地质勘探、自然资源管理、国防、应急响应、电信、传媒、交通、城市规划以及国土安全等领域。
模型概述
Ganos Raster数据模型主要由以下几个元素构成：
Raster：泛指一份栅格数据。例如，一个景遥感影像、一个TIFF文件等。
Tile：数据分块，为一系列像素的集合。Tile为Raster对象在数据库中存储的基本单元，一般每个Tile包含256x256个像素值。
Band：由多个Tile构成的2D栅格数据图层，每个Tile拥有一个行列号。
Cell：表示Tile中的一个像素，可以拥有不同的数据类型。例如Byte，Short，Int，Double等。
Pyramid：通过逐级抽稀的栅格金字塔，方便快速显示。每个Pyramid包含不同的层级，每个层级对应一个Layer，第0层代表原始数据。
Metadata：栅格的存储元数据，包括：空间范围、投影类型、像素类型等。
如上图所示，Ganos Raster采用了一种简单而高效的通用栅格数据模型来管理专题数据和遥感影像数据。一幅栅格数据（Image）在数据库中以栅格对象（Raster）形式进行存储。Raster对象逻辑上由若干可以表示为2D栅格图层的波段（Band）组成，各个Band的信息会在入库过程中从原始影像数据读取。Raster是以数据块（Tile）为基本存储单元进行存储和管理的，Tile的尺寸默认为256x256像素，但也可以由用户进行定义。每个Tile可以包含一个或者多个Band。Tile中的一个像素由一个像素单元（Cell）表示。每一个Raster对象都有对应的元数据（Metadata）信息，如图幅范围（Extent）、数据类型、投影信息、行列号等。 如果对栅格数据创建了金字塔模型，则一个Band会包含多个层级（level）的金字塔数据。
功能优势
相较于开源空间数据库PostGIS自带的Raster插件，Ganos Raster在业务适配、存储成本、计算能力方面都有明显的优势，主要体现在以下几方面：
Ganos Raster存储结构更加贴近业务需求。
与PostGIS Raster完全栅格化的存储方式不同，Ganos Raster提供了面向对象的存储结构，一幅栅格数据（影像，图像、DEM等）导入后直接对应Ganos Raster的一行记录，对用户而言形成了一对一的关系，清晰明了，单个对象的容量没有限制，一行记录可以保存大到1 TB以上的超大影像/栅格。Ganos Raster屏蔽了直接对Tile进行操作的方式，可完整记录栅格对象的元数据信息，并且与时序高度关联，可以与业务进行更好的衔接。
Ganos Raster支持数据降本存储。
由于独特的结构设计，针对海量影像分析这类存储成本较大的场景，Ganos Raster支持将栅格元数据存储在库内，栅格属性数据存储在更为廉价的对象存储上。这种情况下依然支持对栅格进行各类空间分析操作，同时可大幅度降低存储成本。
Ganos Raster拥有更多时空算子。
除了支持传统的栅格空间关系判断、栅格金字塔、栅格像素值、栅格属性、栅格图像处理等操作外，Ganos Raster还支持多种独特的栅格统计、栅格代数操作，同时支持专业的图像匀色算法和海量栅格概视图实现渲染加速。
面雨量分析应用场景
场景描述
面雨量是描述整个区域（流域）内单位面积上平均降水量的物理量，能够较为客观地反映该区域的降水情况。面雨量是洪水预报中极为重要的指标，为各级政府在防汛抗洪及水库调蓄方面提供了依据，同时也为防灾抗灾和经济建设提供参考。由于面雨量作为区域性指标，相关管理部门通常只能获取观测站点提供的点位监测数据，因此需要构建从点位数据到面指标的计算链路，并提供相关附加分析数据以辅助决策。
需求分析
为实现上述场景，需要考虑基于离散监测点进行空间插值，以构建格网数据。随后，基于该网格数据开展等值线/面追踪，同时计算任意矢量面所覆盖范围内的相关统计信息，并求取总体面雨量。
数据导入：使用FDW扩展实现矢量数据的入库，完成矢量观测点数据的入库。
数据插值：使用ST_InterpolateRaster函数基于矢量观测点插值为自定义大小格网，用于存储Ganos Raster。
等值线/面：使用ST_Contour函数基于上述格网数据追踪等值线或面，以分析总体雨量的分布情况。
空间统计：根据不同格网中的雨量分布进行分类，使用ST_Statistics函数统计各类别的最大雨量、最小雨量、中位数及总雨量等信息。
面雨量计算：基于任意矢量边界与格网，计算出矢量所包含的格网。针对穿透的格网，构建相应的统计规则，最终得出该矢量面的总体雨量数据。
技术实现
安装Ganos插件
在目标数据库中安装ganos_raster和ganos_fdw插件。
CREATE EXTENSION ganos_raster CASCADE;
CREATE EXTENSION ganos_fdw CASCADE;导入数据
准备shapefile格式的点要素图层文件point.shp和面要素图层文件polygon.shp，其中每个点要素具有一个id（字符串）属性和一个降雨量pre（浮点型）属性。将图层叠加至QGIS中显示，如下图所示，数字表示每个点的监雨量属性值：
通过Ganos 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)
空间插值
Ganos提供空间插值函数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;
等值线/面追踪
通过插值获得栅格对象后，可以利用Ganos提供的栅格处理函数进行多种处理和分析。以下将以等值线/面为例进行详细介绍。
生成等值线
Ganos的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中叠加配色后，效果如下：
雨量分析统计
以下将展示如何利用前述通过降雨量观测值插值获得的栅格对象，计算任意多边形的面雨量。面雨量是指在特定区域内，降雨总量与该区域面积的比值，通常以毫米/平方米为单位表示。因此，首先计算每个像素覆盖的多边形区域的面积，随后将该像素值乘以相应面积并进行求和，以获得总降雨量。最后，将总降雨量除以多边形的面积即可得到面雨量。
统计任意区域降雨量统计信息。Ganos提供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)
总结
Ganos栅格引擎提供了一整套针对格网数据的导入、存储、分析与可视化支撑能力，帮助高效解决遥感影像管理、DEM数据分析、格网数据分析等业务，同时还具备GPU运算等高阶能力，以支持不同的应用需求。Ganos栅格引擎在面雨量计算场景中的应用，已有效支撑水利部气象降水会商系统的运行，并将该系统的分析计算效率提升了十倍以上。目前，Ganos栅格引擎在水利、自然资源、气象、环境保护、应急等领域具备全面的能力供给与应用案例，能够为多家客户的空天大数据管理类应用提供稳定、高效且健全的时空基底保障。