前言
PolarDB是阿里云自研的云原生关系型数据库，在存储计算分离架构下，利用了软硬件结合的优势，为用户提供具备极致弹性、高性能、海量存储、安全可靠的数据库服务。其中，PolarDB PostgreSQL 版100%兼容 PostgreSQL，高度兼容Oracle语法。Ganos是PolarDB PostgreSQL版提供的新一代云原生时空数据库引擎，在100%兼容PostGIS开源空间数据库基础上，同时具备几何、栅格、轨迹、表面网格、体网格、3D实景、点云、路径、地理网格、快显10大核心引擎完整能力，为数据库构建了面向物理世界时空多模多态数据的混合存储、查询、分析一体化能力。
本文介绍的Ganos栅格数据提取恶劣天气范围，结合Flink进行实时预警能力，依托阿里云云原生关系型数据库PolarDB PostgreSQL版建设输出。
业务场景
远洋船只有很多航线选择，比如最短航线、最经济航线、最安全航线、分段燃油航线等，一旦航线确定下来，船只会沿着既定的航线行驶，且行驶的时间一般都很长。
海洋上气候变化莫测，恶劣天气会对船只造成重大影响，如风浪过大会导致燃油经济效率降低，如遭遇台风、海啸等极端天气甚至会带来财产和生命损失。
因此利用海洋气象预报数据，进行恶劣天气预警，在船只靠近恶劣天气区域时提供信息提示，指导其进行规避，对于航船的安全性，有重要的经济价值和现实意义。
最佳实践
本文以Ganos 从气象数据中提取恶劣天气范围，并实现实时电子围栏能力为例，展示PolarDB在航海恶劣天气预警中的应用。
技术实现
使用PolarDB Ganos 栅格模型对天气预测数据进行管理，提取出恶劣天气空间范围，写入到PolarDB中作为电子围栏表，并结合Ganos的实时电子围栏能力，实现恶劣气象的实时预警能力。
建议配置
为了得到良好的体验，建议使用以下配置：
项目 |
推荐配置 |
PolarDB 版本 |
标准版 兼容PostgreSQL 14 |
CPU |
> 4 Core |
内存 |
>16 GB |
磁盘 |
>100GB (AUTOPL) |
版本 |
>1.1.39 |
Ganos版本 |
>= 6.3 |
气象数据入库
天气预测时会针对某种特定的要素，如风向，风速等进行间隔时间的连续预测，形成一个时间序列，每个时间序列对应到一个时间段，信息会记录到元数据中。
Ganos中每个要素存储为一个栅格对象，每个栅格对象包含若干的波段，每个波段对应到一个时间。
例如，针对一周的每3小时预测数据，会形成 56个时间序列，对应到56个波段。
创建以下表来进行气象数据的存储：
CREATE TABLE weather_table( id serial, name text, -- 文件名称 element text, -- 要素名称 rast raster, -- 栅格对象 ts timestamp, -- 开始时间 te timestamp); -- 结束时间
NetCDF（网络公用数据格式）是一种用来存储温度、湿度、气压、风速和风向等多维科学数据（变量）的文件格式，是一种常见的天气预测数据格式。Ganos可实现对于NetCDF数据格式的直接入库操作（注意： OSS和数据库需要在同一个Region中，并使用内部地址的endpoint）。
普通要素
普通要素如浪高等，入库时不需要进行预处理操作。普通要素可以直接通过ST_ImportFrom 函数进行入库操作，入库过程中不做任何数据修改。
INSERT INTO weather_table(name, element, rast) VALUES ('2023010100_0p25', 'wave_height', ST_SetSRID(ST_ImportFrom('t_chunk_2023010100_0p25','oss://<access_id>:<secrect_key>@[<Endpoint>]/<bucket>/path_to/2023010100_0p25.nc:wave_height'), 4326));
分量要素
对于分量要素，如风向数据，其存的并不是一个方向角度和一个风力大小，而是在南北方向（U）、东西方向（V）两个方向上的风力值，通过分量的方式确定风力和风向，因此在入库时需要进行转换操作。
speed = sqrt (U*U + V*V)
在入库时使用ST_ImportFrom函数先导入到临时表，最终通过ST_MapAlgebra函数进行地图代数运算获得风速信息：
-- create temp table CREATE TEMP TABLE tmp_raster_table(id integer, name text, rast raster); SELECT st_createChunkTable('tmp_chunk_table', true); INSERT INTO tmp_raster_table VALUES (1, '2023010100_0p25', ST_SetSRID(ST_ImportFrom('tmp_chunk_table','oss://<access_id>:<secrect_key>@[<Endpoint>]/<bucket>/path_to/2023010100_0p25.nc:wind_u'), 4326)), (2, '2023010100_0p25', ST_SetSRID(ST_ImportFrom('tmp_chunk_table','oss://<access_id>:<secrect_key>@[<Endpoint>]/<bucket>/path_to/2023010100_0p25.nc:wind_v'), 4326));
-- 此处仅示例6个波段，更多波段可自行添加 WITH foo AS ( SELECT 1 AS rid, rast FROM tmp_raster_table where id = 1 -- u UNION ALL SELECT 2 AS rid, rast FROM tmp_raster_table where id = 2 -- v ) INSERT INTO weather_table(name, element, rast) SELECT '2023010100_0p25', 'wind_speed', ST_SetSRID(ST_MapAlgebra( ARRAY(SELECT rast FROM foo ORDER BY rid), '[{"expr":"sqrt([0,0]**2 + [1, 0]**2)","nodata": true, "nodataValue":-999}, {"expr":"sqrt([0,1]**2 + [1, 1]**2)","nodata": true, "nodataValue":-999}, {"expr":"sqrt([0,2]**2 + [1, 2]**2)","nodata": true, "nodataValue":-999}, {"expr":"sqrt([0,3]**2 + [1, 3]**2)","nodata": true, "nodataValue":-999}, {"expr":"sqrt([0,4]**2 + [1, 4]**2)","nodata": true, "nodataValue":-999}, {"expr":"sqrt([0,5]**2 + [1, 5]**2)","nodata": true, "nodataValue":-999}]', '{"chunktable":"t_chunk_2023010100_0p25"}' ),4326);
计算时间范围
每一个栅格数据都有对应的时间信息，一般以元数据的方式存储在 time#units 和 NETCDF_DIM_time_VALUES 属性中，其中time#units 表示的是开始时间，如 hours since 2023-01-01 00:00:00， NETCDF_DIM_time_VALUES 对应于每个波段的小时数。利用这两个元数据可以计算出对应的开始时间和结束时间：
UPDATE weather_table SET ts = replace(ST_metadata(rast, 'time#units'), 'hours since ', '')::timestamp + (ST_metadata(rast, 'NETCDF_DIM_time_VALUES')::int[])[1] * interval '1 hour', te = replace(ST_metadata(rast, 'time#units'), 'hours since ', '')::timestamp + (ST_metadata(rast, 'NETCDF_DIM_time_VALUES')::int[])[ST_Numbands(rast)] * interval '1 hour' WHERE ts is NULL and te is NULL;
如果我们知道这是一个三小时的预测数据，那可以进一步简化为：
UPDATE weather_table SET ts = replace(ST_metadata(rast, 'time#units'), 'hours since ', '')::timestamp, te = replace(ST_metadata(rast, 'time#units'), 'hours since ', '')::timestamp + interval '3 hour' * （ST_Numbands(rast) + 1） WHERE ts is NULL and te is NULL;
计算恶劣天气范围
恶劣天气的计算符合某些特定要求，比如风速超过一定的阈值或者浪高超过一定的阈值，或者需要具备多种组合条件。
创建函数find_bad_weather_era来获取某种特定要素在特定时间下的区域：
CREATE OR REPLACE FUNCTION find_bad_weather_era(elem text, v float8 ,t timestamp) RETURNS geometry AS $$ DECLARE v_id integer; v_rast raster; v_band integer; v_area geometry; BEGIN -- 获取栅格数据对象 SELECT rast,id into v_rast, v_id FROM weather_table WHERE ts <= t AND te > t and element = elem; RAISE NOTICE 'id = %', v_id; -- 计算波段 SELECT EXTRACT(epoch FROM (t - ts))/3600/3 into v_band FROM weather_table WHERE id = v_id; RAISE NOTICE 'band = %', v_band; -- 计算区域 WITH tmp AS ( SELECT (ST_PixelAsPolygons(v_rast, v_band)).* ) SELECT ST_Union(geom) into v_area FROM tmp WHERE value >= v; return v_area; END; $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT PARALLEL SAFE;
例如， 需要获取时间为 '2023-01-01 3:00:00'， 风力 > 10m/s 的恶劣气象范围， 使用SQL：
SELECT find_bad_weather_era('wind_speed', 10, '2023-01-01 3:00:00');
结果如下图所示
如需要获取海浪 > 4.5 米的恶劣气象，可使用以下SQL：
SELECT find_bad_weather_era('wave_height', 4.5, '2023-01-01 3:00:00');
如果需要满足以上两种与关系，可将两个几何对象进行ST_Intersection 操作；
SELECT ST_INTERSECTION( find_bad_weather_era('wind_speed', 10, '2023-01-01 3:00:00'), find_bad_weather_era('wave_height', 4.5, '2023-01-01 3:00:00') );
如果需要或关系，可进行ST_Union 操作
SELECT ST_UNION( find_bad_weather_era('wind_speed', 10, '2023-01-01 3:00:00'), find_bad_weather_era('wave_height', 4.5, '2023-01-01 3:00:00') );
实时电子围栏
可以根据需求，定期将预测的恶劣天气范围结果保存到PolarDB电子围栏表中，结合Ganos实时电子围栏技术，实现恶劣天气预警。
总结
本文重点介绍了PolarDB Ganos使用栅格模型管理气象数据，使用代数计算对栅格数据进行预处理，并根据条件计算出恶劣天气的空间范围，并结合实时围栏计算，实现恶劣天气的预警预报。 Ganos从栅格数据提供标准时空处理框架，计算效率与综合成本均有大规模改善。未来Ganos还将提供更多高效的面向移动对象的实时计算分析能力，推动相关领域的空间信息应用全面走向“在线化”。
试用体验
欢迎访问PolarDB免费试用页面，选择试用“云原生数据库PolarDB PostgreSQL版”，体验PolarDB Ganos的实时时空计算能力。
