关于Ganos

Ganos是阿里云数据库产品事业部联合阿里巴巴达摩院数据库与存储实验室联合共同研发的新一代云原生位置智能引擎，它将时空数据处理能力融入了云原生关系型数据库PolarDB、云原生多模数据库Lindorm、云原生数据仓库AnalyticDB和云数据库RDS PG等核心产品中。Ganos目前拥有几何、栅格、轨迹、表面网格、体网格、3D实景、点云、路径、地理网格、快显十大核心引擎，为数据库构建了面向新型时空多模多态数据的一体化表达、存储、查询与分析计算能力。

本文介绍的栅格引擎能力，依托阿里云云原生关系型数据库PolarDB、云原生数据仓库AnalyticDB和云数据库RDS PG建设输出。

关于Ganos Raster

栅格数据

栅格数据（Raster Data）是将现实空间划分成规则的网格，每一个网格称为一个单元（像元或像素），并在各单元上赋予相应的属性值来表示实体的一种数据形式。

栅格数据通常有两种类型：专题数据和影像数据。

专题数据：每个栅格像元的值可以是一个测量值或分类值，如高程数据（DEM）、污染物浓度、信号强度、降雨量、土地的权属类型或植被类型等。
影像数据：又称遥感影像或遥感图像（Remote Sensing Image），是通过地面遥感、航空遥感或航天遥感平台拍摄的，记录各种地物电磁波大小的胶片或照片，包含如航空影像和遥感卫星影像等。

每一幅栅格数据均带有时间属性和空间属性，我们称之为时空栅格。时空栅格还强调了栅格数据的时态特性，比如管理时间序列的系列栅格数据。

什么是Ganos Raster

Ganos Raster是Ganos针对栅格数据管理和处理的时空数据引擎（插件形式提供），它允许用户使用阿里云云原生数据库PolarDB、云原生数据仓库AnalyticDB和云数据库RDS存储、索引、查询、分析和传输栅格数据及其相关元数据。栅格数据以数据分块（Tile或Block）的形式存储在上述数据库中支持时间和空间查询。此外，Ganos Raster允许多源栅格数据（如遥感、摄影测量和专题地图）之间的融合与分析以及数据服务发布等功能（如TMS或WMTS等）。Ganos Raster可用于包括气象、环境监测、地质勘探、自然资源管理、国防、应急响应、电信、传媒、交通、城市规划以及国土安全等领域。

Ganos Raster数据模型主要包括以下几个元素构成：

Raster：泛指一份栅格数据，如：一景遥感影像、一个TIFF文件等；
Tile：数据分块，为一系列像素的集合。Tile为Raster对象在数据库中存储的基本单元，每个Tile包含若干Cell（像素或像元），且一般Cell个数为256或512；
Band：由多个Tile构成的2D栅格数据图层，每个Tile拥有一个行列号；
Cell：表示Tile中的一个像素，可以拥有不同的数据类型，如Byte，Short，Int，Double等；
Pyramid：通过逐级抽稀的栅格金字塔，方便快速显示。每个Pyramid包含不同的层级，每个层级对应一个Layer，第0层代表原始数据；
Metadata：栅格的存储元数据，包括：空间范围、投影类型、像素类型等；

Ganos Raster解析

如上图所示，Ganos Raster采用了一种简单而高效的通用栅格数据模型来管理专题数据和遥感影像数据。一幅栅格数据（Image）在数据库中以栅格对象（Raster）形式进行存储。Raster对象逻辑上由若干可以表示为2D栅格图层的波段（Band）组成，各个Band的信息会在入库过程中从原始影像数据读取。Raster是以数据块（Tile）为基本存储单元进行存储和管理的，Tile的尺寸默认为256x256像素，但也可以由用户进行定义。每个Tile可以包含一个或者多个Band。Tile中的一个像素由一个像素单元（Cell）表示。每一个Raster对象都有对应的元数据（Metadata）信息，如图幅范围（Extent）、数据类型、投影信息、行列号等。如果对栅格数据创建了金字塔模型，那么一个Band会包含多个层级（level）的金字塔数据。

Ganos Raster的优势

相较于开源空间数据库PostGIS自带的Raster插件，Ganos Raster在业务适配、存储成本、计算能力方面都有明显的优势，主要体现在：

Ganos Raster存储结构更加贴近业务需求。与PostGIS Raster完全栅格化的存储方式不同，Ganos Raster提供了面向对象的存储结构，一幅栅格数据（影像，图像、DEM等）导入后直接对应Ganos Raster的一行记录，对用户而言形成了一对一的关系，清晰明了，单个对象的容量没有限制，一行记录可以保存大到1TB以上的超大影像/栅格。Ganos Raster屏蔽了直接对Tile进行操作的方式，可完整记录栅格对象的元数据信息，并且与时序高度关联，可以与业务进行更好的衔接；
Ganos Raster支持数据降本存储。由于独特的结构设计，针对海量影像分析这类存储成本较大的场景，Ganos Raster支持将栅格元数据存储在库内，栅格属性数据存储在更为廉价的对象存储上。这种情况下依然支持对栅格进行各类空间分析操作，同时可大幅度降低存储成本；
Ganos Raster拥有更多时空算子。除了支持传统的栅格空间关系判断、栅格金字塔、栅格像素值、栅格属性、栅格图像处理等操作外，Ganos Raster还支持多种独特的栅格统计、栅格代数操作，同时支持专业的图像匀色算法和海量栅格概视图实现渲染加速。

面雨量分析应用场景

场景描述

面雨量是描述整个区域（流域）内单位面积上的平均降水量的物理量，能较客观地反映整个区域的降水情况。面雨量是洪水预报十分重要的指标，能为各级政府防汛抗洪及水库调蓄提供依据，并为防灾抗灾和经济建设服务提供参考。

由于面雨量为区域性指标，而往往相关管理部门可以获取到的为观测站点所提供的点位监测数据，需要构建从点位数据至面指标的计算链路，同时可提供相关附加性分析数据辅助决策。

需求分析

为了实现上述场景，需要考虑基于离散监测点进行空间插值，构建格网数据，然后基于网格数据开展等值线/面追踪，同时计算出任意矢量面所覆盖范围内的相关统计信息并计算总体面雨量。

数据导入：使用Ganos fdw扩展实现矢量数据入库，实现矢量观测点数据入库；
数据插值：使用ST_InterpolateRaster函数将矢量观测点插值为自定义大小格网，存储Ganos raster；
等值线/面：使用ST_Contour函数基于上述格网数据追踪等值线/面，分析总体雨量的分布情况；
空间统计：按不同格网中的雨量分布进行分类，使用ST_Statistics函数根据类别统计出最大雨量、最小雨量、中位数、总雨量等信息；
面雨量计算：基于任意矢量边界与格网计算出矢量包含的格网，对于穿透的格网构建统计规则，最终计算出该矢量面的总体雨量数据。

技术实现

创建Ganos扩展

连接数据库，执行以下SQL语句，创建ganos_raster和fdw扩展。

CREATE EXTENSION ganos_raster CASCADE;
CREATE EXTENSION ganos_fdw CASCADE;

导入数据

准备shapefile格式的点要素图层文件point.shp和面要素图层文件polygon.shp，其中每个点要素拥有一个id（字符串）属性和一个降雨量pre（浮点型）属性。将图层叠加到QGIS显示如下，数字为每个点的监雨量属性值：

用户可以通过Ganos FDW( 参考手册 )将shapefile文件以外表的形式从OSS映射到数据库中，然后进一步通过CREATE方式创建数据库表并插入数据，具体SQL语句如下：

-- 创建fdw服务
-- ak和ak_secret分别为用户OSS服务的AccessKeyId和AccessKeySecret
CREATE SERVER myserver
  FOREIGN DATA WRAPPER ganos_fdw
  OPTIONS (
    datasource 'OSS://ak值:ak_secret值@endpoint名称/bucket名称/子目录',
    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类型的对象。空间插值函数也支持并行操作，通过配置并行度提高执行效率。

第1步：创建表格raster_table，其中包一个raster类型的字段，用于保存插值后的栅格数据。

CREATE TABLE IF NOT EXISTS raster_table
  id integer, 
  rast raster -- raster数据类型，用于保存插值后的栅格对象
);

第2步：调用ST_InterpolateRaster函数对point_table表中空间点数据通过反距离加权(IDW)方法进行插值，并将结果插入到raster_table中。注意这里我们使用ST_MakePoint方法创建了一个带属性值（降雨量）的Point对象作为插值的输入参数：

INSERT INTO raster_table(id, rast) VALUES(
  (SELECT ST_InterpolateRaster(
    ST_Collect(ST_MakePoint(ST_X(geom),ST_Y(geom),pre)),
    '{"method":"IDW","radius":2.0,"max_points":"30","min_points":"1"}',
    '{"chunktable":"rbt","celltype":"32bf"}')
  FROM point_table));

当插值点数量比较大时候，用户可以在执行插值函数执行前开启并行来实现函数加速：

-- 配置并行度为4
SET ganos.parallel.degree = 4;

第3步: 用户可以通过ST_ExportTo将生成的raster对象导出为文件形式方便查看，这里我们以COG格式导出到OSS：

SELECT ST_ExportTo(rast, 
  'COG',
  'OSS://ak值:ak_secret值@endpoint名称/bucket名称/子目录/idw.tif') 
  FROM raster_table WHERE id=1;

将文件文件下载并加载到QGIS效果如下所示：

第4步: 调用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://ak值:ak_secret值@endpoint名称/bucket名称/子目录/idw_clip.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属性并指定为ture，则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中叠加配色后效果如下：

雨量分析统计

下面我们来展示如何利用前面通过降雨量观测值插值获取的栅格对象来计算任意多边形的面雨量。面雨量是指在一个区域内,降雨的总量与该区域面积的比值,通常以毫米/平方米为单位表示，所以本案例我们需要先计算出每个像素覆盖的多边形区域的面积，然后乘以该像素值并求和，从而得到总降雨量，然后除以多边形的面积即可获得面雨量。

第1步：统计任意区域降雨量统计信息

Ganos提供用于栅格数据统计的函数接口ST_Statistics( 参考手册 )，允许用对任意区域(可以是点、线、面等)内的栅格像素值进行统计。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个等。

第2步：裁剪栅格对象

参考前面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;

为了获取最准确的面雨量统计信息，我们需要遍历该栅格数据的每一个像素，并获取该像素覆盖区域的矩形要素和对应的像素值。所以第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

第3步：获取每个像素的面积与降雨量

获取栅格长宽属性后，我们就可以对像素进行遍历，并通过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);
--遍历裁剪栅格数据的每一个像素，获取其覆盖的矩形区域和对应的像素值
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;
$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)

第4步：计算面雨量

最后我们就可以计算面雨量了。考虑到面雨量单位是按照每平方米进行计算，所以在计算每个像素多边形面积之前，需要先通过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)"

基于Ganos的栅格引擎开展区域面雨量分析