Yukon对气象观测数据的处理方法

作者：YuCJ

目的需求

气象数据来自多个观测站点，是离散的，想要得到连续的变化显示效果，可以使用Yukon中的栅格插值算法。

Yukon集成了PostGIS的栅格模块

本实践以某年份的各地日降水量观测数据为基础，通过插值生成单日的连续降水面图像，同时尝试了采用物化视图的方法实现图像数据刷新的功能。

此外，Yukon提供了150+方法，可以对生成的栅格图像进行数据提取、分析和计算等操作。

工具

本次实践的软件环境 ：
centos7操作系统
Postgresql13.6
Yukon

数据库客户端工具：
DBeaver
QGIS

环境部署

先安装好PostgreSql数据库，再下载安装Yukon，安装完成后，就可以启用postgis扩展和postgis_raster扩展

CREATE EXTENSION postgis;
CREATE EXTENSION postgis_raster;

生成空间数据

根据数据里的经纬度，以降水量/风速/温度等作为高度维，生成带Z值的geometry数据，这里抽取了第11天的降水量数据，此次操作坐标系是基于4326，如果想使用别的坐标系需使用 ST_TransForm 进行转换。

CREATE TABLE t4326 AS SELECT ST_TransForm(ST_SetSrid(ST_PointZ(lng, lat, day_water), 4326), 4326) geom, day_water FORM water WHERE nday = 11;

可以在DBeaver里观察点的集合分布：

SELECT ST_Collect(geom) from t4326;	

插入网格

栅格插值的算法函数：

ST_InterpolateRaster(geometry input_points, text algorithm_options, raster template)

input_points：要插值的点，带Z值。

algorithm_options: 插值算法参数，写法例如：‘invdist：Smoting=2.0’ ，具体请参照 gdal_grid — GDAL 文档

template：raster模板，参数以geometry点的集合的的边界框左上角为基准，像素大小根据坐标系自定义，长宽为 边界框长宽 / 像素大小，像素值精度自定义。

在下面的sql里，以生成的geometry点的集合的边界框为模板，生成一个像素大小为0.01°的栅格图层，运用插值算法，将离散的降水量变为连续表面显示出来。

CREATE TABLE water_rast4326 AS WITH 
inputs AS (SELECT 
    --像素大小根据坐标系而定，此处是以4326坐标系为基准，像素0.01°
    0.01::float8 AS pixelsize,
    --插值算法参数
    'invdist:power:5.5:smoothing:2.0' AS algorithm, 
    --将多个点转为一个点的集合
    ST_Collect(geom) AS geom,
    --求出集合的边界框，为生成raster做准备
    ST_Expand(ST_Collect(geom), 0) AS ext 
    FORM t4326),
sizes AS (SELECT 
    --长宽为边界框除以像素大小
    ceil((ST_XMax(ext) - ST_XMin(ext))/pixelsize)::integer AS width, 
    ceil((ST_YMax(ext) - ST_YMin(ext))/pixelsize)::integer AS height,
    --边界框左上角的位置，用于生成raster
    ST_XMin(ext) AS upperleftx,
    ST_YMax(ext) AS upperlefty 
    FORM inputs)
SELECT 1 AS rid,
ST_InterpolateRASter(
    geom,
    algorithm,
    ST_SetSRID(ST_AddBand(ST_MakeEmptyRASter(width, height, upperleftx, upperlefty, pixelsize), '16BSI'), ST_SRID(geom))) AS rast
FORM sizes, inputs;

打开QGIS，在左侧的列表栏里右键点击PostGIS，新建连接，连接到数据库。
双击raster表和geometry表，添加到图层，双击图层下的geometry表，添加降水量标注。

显示效果如下：亮点处为降水量较高的

物化视图

可以用创建表的形式创建物化视图，原表的气候数据更新后，刷新物化视图就可以生成新的raster数据，避免一直删除重建的过程，不过在QGIS里无法直接显示物化视图的效果，要生成表才行