基于PolarDB Ganos的气象数据处理：恶劣气象预警

摘要

本文介绍了在航海恶劣天气预警的场景中，PolarDB Ganos利用栅格模型对基于NetCDF的气象预测数据入库，利用数据库函数计算恶劣天气范围，并将其存储在PolarDB中。通过Ganos实时电子围栏功能，对船只位置并与恶劣天气范围对比，实现航海恶劣气象预警功能。

前言

PolarDB是阿里云自研的云原生关系型数据库，在存储计算分离架构下，利用了软硬件结合的优势，为用户提供具备极致弹性、高性能、海量存储、安全可靠的数据库服务。其中，PolarDB PostgreSQL 版100%兼容 PostgreSQL，高度兼容Oracle语法。Ganos是PolarDB PostgreSQL版提供的新一代云原生时空数据库引擎，在100%兼容PostGIS开源空间数据库基础上，同时具备几何、栅格、轨迹、表面网格、体网格、3D实景、点云、路径、地理网格、快显10大核心引擎完整能力，为数据库构建了面向物理世界时空多模多态数据的混合存储、查询、分析一体化能力。

本文介绍的Ganos栅格数据提取恶劣天气范围，结合Flink进行实时预警能力，依托阿里云云原生关系型数据库PolarDB PostgreSQL版建设输出。

业务场景

远洋船只有很多航线选择，比如最短航线、最经济航线、最安全航线、分段燃油航线等，一旦航线确定下来，船只会沿着既定的航线行驶，且行驶的时间一般都很长。

海洋上气候变化莫测，恶劣天气会对船只造成重大影响，如风浪过大会导致燃油经济效率降低，如遭遇台风、海啸等极端天气甚至会带来财产和生命损失。

因此利用海洋气象预报数据，进行恶劣天气预警，在船只靠近恶劣天气区域时提供信息提示，指导其进行规避，对于航船的安全性，有重要的经济价值和现实意义。

最佳实践

本文以Ganos 从气象数据中提取恶劣天气范围，并实现实时电子围栏能力为例，展示PolarDB在航海恶劣天气预警中的应用。

技术实现

使用PolarDB Ganos 栅格模型对天气预测数据进行管理，提取出恶劣天气空间范围，写入到PolarDB中作为电子围栏表，并结合Ganos的实时电子围栏能力，实现恶劣气象的实时预警能力。

建议配置

为了得到良好的体验，建议使用以下配置：

气象数据入库

天气预测时会针对某种特定的要素，如风向，风速等进行间隔时间的连续预测，形成一个时间序列，每个时间序列对应到一个时间段，信息会记录到元数据中。

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 (
SELECT1AS rid, rast FROM tmp_raster_table whereid = 1-- u
UNIONALL
SELECT2AS rid, rast FROM tmp_raster_table whereid = 2-- v
)
INSERTINTO weather_table(name, element, rast)
SELECT'2023010100_0p25', 'wind_speed', 
  ST_SetSRID(ST_MapAlgebra(
    ARRAY(SELECT rast FROM foo ORDERBY 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 ORREPLACEFUNCTION 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,idinto v_rast,  v_id
    FROM weather_table
    WHERE ts <= t AND te > t andelement = elem;

    RAISE NOTICE 'id = %', v_id;

    -- 计算波段
    SELECTEXTRACT(epoch FROM (t - ts))/3600/3into v_band
    FROM weather_table
    WHEREid = 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 
    WHEREvalue >= 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的实时时空计算能力。