很多出国深造的同学，都对国外高校中的计算机教学、使用记忆犹新。国内一般院校的老师很多都是从微软的DOS起步开始捣鼓微型计算机的，基本上对unix系统用的不多。对命令行操作，也停留在dos命令的概念上。最近，一位同学毕业设计遇到了读取天气预报数据并显示在地图上的问题，来请教我，我们一起在linux下摸索了很久，终于搞定了。在学习过程中，参照了几篇前人的文章，帮助很大。

基于GFS数据开发行业气象信息API（I）

基于GFS数据开发行业气象信息API（II）

C++中调用cmd命令行运行脚本处理GDAS数据.

感谢这些朋友的分享！

1 问题背景

该学生毕业设计是做一个天气图层，希望GIS显示出各个区域的温度高低。搜索了很多现有的接口，都不满意，最终在CSDN上“基于GFS数据开发行业气象信息API（I）”，找到了一个免费的高大上数据源，在这个网站，每6个小时会更新一次全球的天气预报，而且是免费的。

GFS应该是全球预报系统的缩写，“GDAS”也是一种天气数据缩写。但这两个东西，距离GIS图层可视化还差了很远很远，实质上，网站发布的是每6小时发布一次的数据文件。这个数据文件存储了很多专题的数据，高度压缩到一个文件里面去。

学生要求我帮助处理的文件名类似gdas.t00z.pgrb2.0p25.f001，“t00”大概表示格林威治时刻0时发布的文件。0p25指的是精度，每0.25度（经纬度）一个点。f001表示的是预测第一个小时的天气。当然，不同精度、不同发布时刻和预测时刻，对应的文件名不同。

该同学的任务是帮助毕业设计小组，提取出一个区域的温度数据，并入库。

2 使用wgrib2读取数据

根据网站的说明，知道需要用一个叫做“wgrib2”的工具进行操作。也是根据官网的指南，发现这个工具可以直接从linux中构建，也可以在windows下安装OpenGrADS，里面直接附带了wgrib2.

2.1 查看数据成分

我们尝试下载一个文件，运行命令行：

[email protected]$wgrib2 gdas.t00z.pgrb2.p25.f000

可以看到数据的所有专题条目，每一行对应了一种数据：

::d=:UGRD:planetary boundary layer:anl:
::d=:VGRD:planetary boundary layer:anl:
::d=:VRATE:planetary boundary layer:anl:
::d=:GUST:surface:anl:
::d=:HGT: mb:anl:
::d=:TMP: mb:anl:
::d=:RH: mb:anl:
::d=:UGRD: mb:anl:
::d=:VGRD: mb:anl:
::d=:O3MR: mb:anl:
::d=:HGT: mb:anl:
::d=:TMP: mb:anl:
……
::d=:PWAT:entire atmosphere (considered as a single layer):anl:
::d=:CWAT:entire atmosphere (considered as a single layer):anl:
::d=:RH:entire atmosphere (considered as a single layer):anl:
::d=:TOZNE:entire atmosphere (considered as a single layer):anl:
::d=:HLCY:- m above ground:anl:
::d=:USTM:- m above ground:anl:
::d=:VSTM:- m above ground:anl:
::d=:PRES:tropopause:anl:
::d=:ICAHT:tropopause:anl:
……
::d=:VWSH:PV=-e- (Km^/kg/s) surface:anl:
::d=:PRMSL:mean sea level:anl:
::d=:WAVH: mb:anl:

包罗万象。可是，这些数据是什么意思呢？查看命令行选项，才发现这个wgrib2真的好好多的参数哦！最终发现了-V和-v参数，可以输出更多的内容。

[email protected]$wgrib2 gdas.t00z.pgrb2.p25.f000 -V

输出：

::vt=:planetary boundary layer:anl:UGRD U-Component of Wind [m/s]:
    ndata=:undef=:mean=.:min=-:max=
    grid_template=:winds(N/S):
    lat-lon grid:( x ) units e- input WE:NS output WE:SN res 
    lat  to - by .
    lon . to  by . #points=1038240

::vt=:planetary boundary layer:anl:VGRD V-Component of Wind [m/s]:
    ndata=:undef=:mean=.:min=-:max=
    grid_template=:winds(N/S):
    lat-lon grid:( x ) units e- input WE:NS output WE:SN res 
    lat  to - by .
    lon . to  by . #points=1038240

……
::vt=: mb:anl:TMP Temperature [K]:
    ndata=:undef=:mean=:min=:max=
    grid_template=:winds(N/S):
    lat-lon grid:( x ) units e- input WE:NS output WE:SN res 
    lat  to - by .
    lon . to  by . #points=1038240

……
::vt=:mean sea level:anl:PRMSL Pressure Reduced to MSL [Pa]:
    ndata=:undef=:mean=:min=:max=
    grid_template=:winds(N/S):
    lat-lon grid:( x ) units e- input WE:NS output WE:SN res 
    lat  to - by .
    lon . to  by . #points=1038240

::vt=: mb:anl:WAVH -Wave Geopotential Height [gpm]:
    ndata=:undef=:mean=:min=:max=
    grid_template=:winds(N/S):
    lat-lon grid:( x ) units e- input WE:NS output WE:SN res 
    lat  to - by .
    lon . to  by . #points=1038240

::vt=:surface:anl:LANDN Land-sea coverage (nearest neighbor) [land=,sea=] [-]:
    ndata=:undef=:mean=.:min=:max=
    grid_template=:winds(N/S):
    lat-lon grid:( x ) units e- input WE:NS output WE:SN res 
    lat  to - by .
    lon . to  by . #points=1038240

不但包含了我们需要的温度参数，还有好多好多东西，比如气压、风力、风向。

2.2 尝试导出实际数据

仔细阅读了命令行说明，发现这个工具竟然可以直接导出为bin或者csv、text

不妨试试看：

$wgrib2 gdas.t00z.pgrb2.p25.f000 -csv .csv -i

而后，直接输入某一行数据，回车，即可生成一个csv文件了。

比如：

$wgrib2 gdas.t00z.pgrb2.p25.f000 -csv .csv -i
::d=:TMP: mb:anl:[回车]

产生了一个csv文件，用wps直接打开。

产生时刻	预报时刻	内容	气压？	纬度（度）	经度（度）	温度（K）
……
2018-05-26 00:00:00	2018-05-26 00:00:00	TMP	400 mb	-0.5	-69	222.6
2018-05-26 00:00:00	2018-05-26 00:00:00	TMP	400 mb	-0.25	-69	222.6
2018-05-26 00:00:00	2018-05-26 00:00:00	TMP	400 mb	-68.75	222.6
2018-05-26 00:00:00	2018-05-26 00:00:00	TMP	400 mb	0.25	-68.75	222.5
……

可见，如果有一个数据库，可以直接录入。不过，由于存在很多组温度数据，对应了不同的大气层压力，也就是海拔，我们如果每一个参数都需要手工的导出，未免太麻烦了。实际上，官方文档中，有例子使用管道直接导出。

3 使用 grep 配合管道批量筛选数据

unix 下的很多程序，都是使用管道进行通信的。使用管道进行通信，可以级联使用，即使用“|”符号前后贯通各个进程。当使用管道时，“|”符号左侧的进程标准输出(stdout)会和符号右侧进程的标准输入（“stdin”）直接对接。这就相当于是把屏幕输出直接送给了下一个进程的输入。

3.1 管道操作

按照官网的例子，我们可以使用：

$wgrib2 gdas.t00z.pgrb2.p25.f000 | grep ':TMP:' |wgrib2 gdas.t00z.pgrb2.p25.f000 -lola :: :: out.csv spread  -i

直接获取东经100~125度（100,100+100×0.25），北纬 20~45度（20,20+100×0.25）的所有温度数据，存储为csv格式。筛选得到感兴趣的数据，有了结构化的数据，故而直接创建表格即可。表格包括时刻、经度、纬度、气压（海拔）、温度几个参数。

3.2 详细解释

这里需要解释的就是linux的管道啦！其实，这个命令行共运行了3个主要进程。

进程1

wgrib2 gdas.t00z.pgrb2.p25.f000

会列出本文件中包含的所有信息，正如2.1里面列出的一模一样。

进程2

会把进程1的输出，作为输入，并筛选出含有温度标记“TMP：”的行。

我们运行

可以看到

::d=:TMP: mb:anl:
::d=:TMP: mb:anl:
::d=:TMP: mb:anl:
……
::d=:TMP:PV=e- (Km^/kg/s) surface:anl:
::d=:TMP:PV=-e- (Km^/kg/s) surface:anl:

进程3

wgrib2 gdas.t00z.pgrb2.p25.f000 -lola :: :: out.csv spread  -i

即执行了多次2.2节的操作，实质上，进程2的输出作为进程3的输入，等于手工敲了好多行。

此外，一些额外的参数用来规定格式、区域：

-lola 参数

规定区域，后面跟着四部分：

经度范围	纬度范围	文件名	格式
开始经度：采集个数：采集间隔	开始纬度：采集个数：采集间隔	输出文件	文件格式
100:100:0.25	20:100:0.25	out.csv	spread

-i 参数

表示从stdin输入要导出的内容。

导出结果

可以看到，导出的数据如下：

longitude	latitude	value
100.000000	20.000000	259.5
100.250000	20.000000	259.6
……
124.500000	44.750000	264.9
124.750000	44.750000	264.9

文本文件每组一万行，即100×100个采集经纬度位置。共有好多组，分别对应了筛选出的各组数据。我们把这个交给毕业设计小组的组员，就能够导入PostgreSQL数据库中啦。

4 使用octave 可视化温度数据

目前，学生的GIS还没做好。我建议他使用Qt+fcgi的方式直接生成半透明的图片叠加到OpenStreetMap上去。现在，我们先使用octave对温度进行读取，可以直接绘制出全球的温度等值线。

为了方便octave读取，我们输出为二进制格式，还是用管道一步到位。

$wgrib2 gdas.t00z.pgrb2.p25.f000 | grep ':TMP:' |wgrib2 gdas.t00z.pgrb2.p25.f000 -lola :: :: out.bin bin  -i

这样会生成一个含有所有温度信息的二进制文件。二进制文件的格式是

字节长度1 ，内容，字节长度1，字节长度2，内容，字节长度2，……

内容直接就是float32类型的矩阵，由于我们导出的是一个区域，实际上有100x100个点。绘制脚本仅仅导出一部分数据进行绘制，否则太多了。

clear
lats = ::;
lons = ::;
[mx,my] = meshgrid(lons,lats);

fid = fopen('out.bin','rb');
asize = fread(fid,[],'uint32');
allv = fread(fid,[,],'float32');
ct = ;
while (asize==)
  ct = ct + ;
  if (ct>= && ct <=)
    subplot(,,ct-);
    mesh(mx,my,allv);  
  endif
  asize = fread(fid,[],'uint32');
  asize = fread(fid,[],'uint32');
  allv = fread(fid,[,],'float32');
endwhile
fclose(fid);

输出图像是几个海拔高度下的温度分布，单位为开尔文（K）：

管道操作——为GIS准备GDAS气象数据1 问题背景2 使用wgrib2读取数据3 使用 grep 配合管道批量筛选数据4 使用octave 可视化温度数据4 后记

4 后记

在帮助学生完成了后续的可视化工作后，还是觉得意犹未尽。不知道学习气象的同学知不知道国内气象数据的开发情况，特别是有没有这样的开放数据网站。学生给我看的国内的天气预报接口都太花哨、却不够深入、方便——做地理信息的展示，原始数据是最合适的选择。

4.1关于数据网站

看到这个网站的美工，不由得想起90年代的门户。虽然很土，但是配合了如此的内容，真的有一种老派JazzBand的感觉。如果稍微深入阅读一下网站的说明，就会发现这个数据源背后，是一群从80年代就开始，用fortran写代码、在unix系统里做研究的超级神人，以及他们带出的学生。能够把如此多的科学参数，整合到一起，有序的管理起来，同时，开发出非常灵活、有效的工具链进行发布，本身就非常厉害了。更可贵的是，文档也非常通俗易懂，这种科学精神真的值得我们学习。

4.2留学生的计算机技能准备

出国留学的同学，一定会发现国外大学与国内的不同。有些大学论文要求用latex，周围同学做笔记，用的也是libreoffice等等开源文本处理器。计算实验用Matlab学生版是有的，也有VC++ 可用，但学习一下python, Octave，GNU C, Linux绝对没有坏处——这样想抄别人的也好抄。

管道操作——为GIS准备GDAS气象数据1 问题背景2 使用wgrib2读取数据3 使用 grep 配合管道批量筛选数据4 使用octave 可视化温度数据4 后记

1 问题背景

2 使用wgrib2读取数据

2.1 查看数据成分

2.2 尝试导出实际数据

3 使用 grep 配合管道批量筛选数据

3.1 管道操作

3.2 详细解释

进程1

进程2

进程3

-lola 参数

-i 参数

导出结果

4 使用octave 可视化温度数据

4 后记

4.1关于数据网站

4.2留学生的计算机技能准备

继续阅读

Apache (You don't have permission to access / on this server.）

debian9升级4.9.0内核到4.19.2内核过程

centOS7 配置 vsftpd 虚拟用户及权限Vsftpd配置虚拟用户及权限

linux-svn卸载与安装

vsftp虚拟多用户多权限一键部署脚本

Ubuntu14.04 LTS下安装mongodb

Nginx服务优化（1）——隐藏版本号、修改用户与组、网页缓存时间、日志切割、连接超时一、隐藏版本号二、修改用户与组三、配置Nginx网页缓存时间四、实现Nginx日志分割五、配置Nginx实现连接超时六、补充关于时间日期的命令

httpd服务的部署、启动、配置和简单优化一、部署二、启动三、配置文件

配置网页内容访问

手动安装Intel network I217-LM网卡的Linux驱动

禁止ubuntu系统弹出报错界面

Ubuntu Linux下Apache的配置文件

samba服务器的功能

【Linux】UDP广播报文接收速率问题

Linux设备模型（中）之上层容器

PowerPC平台 Linux移植三