格点数据应用之一 数据存储篇（PG）

格点数据是一种常见的数据类型，典型如气象领域的NetCDF与Grib/Grib2文件，但是由于格点数据比较冷门，专业搞数据库设计的技术人员搞不懂也不知道如何设计存储，而另一方面搞气象应用开发的技术人员大部分情况下不懂数据库还在使用格点文件进行业务应用开发，因此特意写下本篇文章，能对专业搞数据库的和专业搞气象应用的技术人员有所启发。

本文基于PostgreSQL数据库进行设计，PG数据库是一个TP与AP都比较强势的混合型数据库，特性多功能强黑科技丰富。本文主要技术思想参考自《基于云上分布式NoSQL的海量气象数据存储和查询方案》，并结合PG实际情况进行稍微改良。

格点数据是典型的多维数据集，打开任何格点数据文件，如下图示意：

netcdf（nc）文件文件都会显示列出（要素），时间，空间维度信息，不同维度组合能得到不同的格网平面数据，因此一个格点文件是由多个不同格网平面数据组成的。概括来说，格点文件维度信息主要由“四维”构成，说明如下：

结论：格点文件是多维数据集，文件内部根据维度组合由多个格网平面数据组成，通过“四维”能唯一描述某一个格网平面。

首先贴个格网示意图：

格网本身是有两部分组成：元数据+格点数组，元数据会定义格网的维度特征，例如包括格网的起点经纬度Start，终点经纬度End，除此之外包括一些其他重要信息，如每个格点（小方格就是格点）的宽高（一般都相等）res_x，res_y，宽就是水平方向的分辨率，高就是垂直方向的分辨率，一个格点文件内的格网平面会共享部分重要元数据，列表信息如下：

同时，不同格网平面自身会根据“四维”参数不同，有自己的独享元数据，列表信息如下：

元数据与格点值数组之间存在很多计算关系。

观察上图可知，格网的宽高是由地理范围和格点分辨率决定的，计算公式如下：

格网的行列号与格点数组索引序号存在如下计算关系：

一个格网平面是由很多格点单元（小正方形）组成的，这种小单元拥有地理位置属性，即输入任意经纬度，都能通过元数据，计算出其对应格网行列号，再通过上一个公式获取格点值，具体计算关系如下：

格点数据在业务中应用场景主要分三部分：

很多生活类网站诸如“未来7天天气”这种功能，如下图：

前文可知，格点文件是由多个格网平面组成，格网平面由多个格点单元组成，格点单元存在地理属性，所谓穿透查询，就是通过输入任意经纬度，根据元数据转换成格网行列号，然后获取相关格网平面在该行列号处的格点值。细分有：1 根据预报时刻的穿透查询，输入的是当前所在城市的经纬度，输出连续若干天的格点值对应的天气情况；2 根据起报时间的历史天气查询，如“历史上的今天的天气怎么样？”； 3 根据高度不同，查询从高空到低空的对应地理位置天气情况等。

属于可视化范畴，具体参考《格点数据应用之二 WebGIS等值面可视化篇》。

如国家电网需要根据暴雨等级等值面分析受影响的高压杆塔等，属于GIS空间分析范畴，这一部分是行业中最有用的气象+行业的一种应用场景，具体可参考本公众号Spatial Data后续文章《格点数据应用之三 空间分析性能优化篇》。

本篇主要阐述格点穿透查询，由上文阐述可知，穿透查询由于查询的条件不同，可能涉及多个气象格点文件，通常每个格点文件都比较大，如果使用格点文件进行业务应用时，通常需要先打开若干业务相关格点文件，把较大数据加载到内存取数，IO与CPU开销都很大，这种做法是很不合理的，因此本文结合PG做了存储设计与性能测试。

前文所述，一个格点文件，通常包括很多个不同时间不同高度的格网，而一个格网，是由元数据+格点数组组成的。采用连续数组是一种空间连续数据的一种压缩形式，避免反复描述数据位置时间等增大文件体积，所以在设计上我们仍将保留这种格点文件设计形式。

存储设计上将一份格网设计成三张表：数据集表（dataset)，格网表（grid），数据表（data），设想一个格点文件是一个数据集，数据集表主要描述通用元数据，格网表主要描述“四维”独享元数据，数据表采用PG数组类型存储格点数组值。各个表之间关系如下：一个数据集对应多个格网，每个格网都会标记所属的数据集id，一个格网对应数据表里多个格点数组，这里说明下，格网文件是用一个数组存的数据，而数据库将这个大的数组拆分成连续的小数组存储，避免单列太到影响性能。

这里详细说明下为何拆分数据数组，由于穿透查询只针对局部某个格点，从小数组中获取数值比从大数组中获取数值能极大降低IO开销。本文实际测试通常拆分小数组以2000到2500长度为宜，该列大小小于8kb，PG库Page默认8kb，如果数据超过8kb也会增加扫描开销，那么，假设100*100的格网，共计10000个格点值，这里如果按照2500一组拆分，实际数据表就会存四行数据，每个记录所属格点id，分组序号，分组数组。

详细表结构设计如下：

那么穿透查询，输入任意经纬度，需要先转换行列号，再根据行列号和格网分块数量，计算对应的chuck_idx值和datas数组的值即可，计算公式如下：

示例：查询某一数据集id为xxx，物理量为x，qbsj为n,ybsk为[xx1,xx2,...xxn]，高度默认0的数据天气情况，伪代码分两步：第一步：根据数据集id获取元数据，根据输入经纬度计算chuck_idx与datas_idx

第二部分，根据chuck_idx与datas_idx去数据库捞数据：

通常将第一部分第二部分合并业务逻辑写成PG的一个function，对外提供业务服务。笔者实际测试情况是，格点穿透查询基本在1-2ms之间，速度非常快。