很多出国深造的同学,都对国外高校中的计算机教学、使用记忆犹新。国内一般院校的老师很多都是从微软的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 查看数据成分
我们尝试下载一个文件,运行命令行:
user@localhost$wgrib2 gdas.t00z.pgrb2.0p25.f000
可以看到数据的所有专题条目,每一行对应了一种数据:
1:0:d=2018052600:UGRD:planetary boundary layer:anl:
2:511774:d=2018052600:VGRD:planetary boundary layer:anl:
3:1012368:d=2018052600:VRATE:planetary boundary layer:anl:
4:1537417:d=2018052600:GUST:surface:anl:
5:2075143:d=2018052600:HGT:1 mb:anl:
6:2548248:d=2018052600:TMP:1 mb:anl:
7:2771646:d=2018052600:RH:1 mb:anl:
8:2895954:d=2018052600:UGRD:1 mb:anl:
9:3511519:d=2018052600:VGRD:1 mb:anl:
10:3847835:d=2018052600:O3MR:1 mb:anl:
11:4250299:d=2018052600:HGT:2 mb:anl:
12:4722967:d=2018052600:TMP:2 mb:anl:
……
274:144801219:d=2018052600:PWAT:entire atmosphere (considered as a single layer):anl:
275:145242706:d=2018052600:CWAT:entire atmosphere (considered as a single layer):anl:
276:145565348:d=2018052600:RH:entire atmosphere (considered as a single layer):anl:
277:145840237:d=2018052600:TOZNE:entire atmosphere (considered as a single layer):anl:
278:146289807:d=2018052600:HLCY:3000-0 m above ground:anl:
279:146836402:d=2018052600:USTM:6000-0 m above ground:anl:
280:147639100:d=2018052600:VSTM:6000-0 m above ground:anl:
281:148427006:d=2018052600:PRES:tropopause:anl:
282:149355243:d=2018052600:ICAHT:tropopause:anl:
……
350:189417192:d=2018052600:VWSH:PV=-2e-06 (Km^2/kg/s) surface:anl:
351:189746442:d=2018052600:PRMSL:mean sea level:anl:
352:190290484:d=2018052600:5WAVH:500 mb:anl:
包罗万象。可是,这些数据是什么意思呢?查看命令行选项,才发现这个wgrib2真的好好多的参数哦!最终发现了-V和-v参数,可以输出更多的内容。
user@localhost$wgrib2 gdas.t00z.pgrb2.0p25.f000 -V
输出:
1:0:vt=2018052600:planetary boundary layer:anl:UGRD U-Component of Wind [m/s]:
ndata=1038240:undef=0:mean=0.33333:min=-30.1603:max=32.6397
grid_template=0:winds(N/S):
lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48
lat 90.000000 to -90.000000 by 0.250000
lon 0.000000 to 359.750000 by 0.250000 #points=1038240
2:511774:vt=2018052600:planetary boundary layer:anl:VGRD V-Component of Wind [m/s]:
ndata=1038240:undef=0:mean=0.495147:min=-36.3917:max=38.5083
grid_template=0:winds(N/S):
lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48
lat 90.000000 to -90.000000 by 0.250000
lon 0.000000 to 359.750000 by 0.250000 #points=1038240
……
6:2548248:vt=2018052600:1 mb:anl:TMP Temperature [K]:
ndata=1038240:undef=0:mean=257.683:min=228.2:max=278.5
grid_template=0:winds(N/S):
lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48
lat 90.000000 to -90.000000 by 0.250000
lon 0.000000 to 359.750000 by 0.250000 #points=1038240
……
351:189746442:vt=2018052600:mean sea level:anl:PRMSL Pressure Reduced to MSL [Pa]:
ndata=1038240:undef=0:mean=101040:min=95558.1:max=104411
grid_template=0:winds(N/S):
lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48
lat 90.000000 to -90.000000 by 0.250000
lon 0.000000 to 359.750000 by 0.250000 #points=1038240
352:190290484:vt=2018052600:500 mb:anl:5WAVH 5-Wave Geopotential Height [gpm]:
ndata=1038240:undef=0:mean=5532.88:min=4737.75:max=5940.99
grid_template=0:winds(N/S):
lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48
lat 90.000000 to -90.000000 by 0.250000
lon 0.000000 to 359.750000 by 0.250000 #points=1038240
353:190617698:vt=2018052600:surface:anl:LANDN Land-sea coverage (nearest neighbor) [land=1,sea=0] [-]:
ndata=1038240:undef=0:mean=0.338136:min=0:max=1
grid_template=0:winds(N/S):
lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48
lat 90.000000 to -90.000000 by 0.250000
lon 0.000000 to 359.750000 by 0.250000 #points=1038240
不但包含了我们需要的温度参数,还有好多好多东西,比如气压、风力、风向。
2.2 尝试导出实际数据
仔细阅读了命令行说明,发现这个工具竟然可以直接导出为bin或者csv、text
不妨试试看:
$wgrib2 gdas.t00z.pgrb2.0p25.f000 -csv 1.csv -i
而后,直接输入某一行数据,回车,即可生成一个csv文件了。
比如:
$wgrib2 gdas.t00z.pgrb2.0p25.f000 -csv 1.csv -i
125:55719331:d=2018052600:TMP:400 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 | 0 | -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.0p25.f000 | grep ':TMP:' |wgrib2 gdas.t00z.pgrb2.0p25.f000 -lola 100:100:0.25 20:100:0.25 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.0p25.f000
会列出本文件中包含的所有信息,正如2.1里面列出的一模一样。
进程2
grep ':TMP:'
会把进程1的输出,作为输入,并筛选出含有温度标记“TMP:”的行。
我们运行
wgrib2 gdas.t00z.pgrb2.0p25.f000 | grep ':TMP:'
可以看到
6:2548248:d=2018052600:TMP:1 mb:anl:
12:4722967:d=2018052600:TMP:2 mb:anl:
18:6861957:d=2018052600:TMP:3 mb:anl:
……
341:185271533:d=2018052600:TMP:PV=2e-06 (Km^2/kg/s) surface:anl:
347:187887025:d=2018052600:TMP:PV=-2e-06 (Km^2/kg/s) surface:anl:
进程3
wgrib2 gdas.t00z.pgrb2.0p25.f000 -lola 100:100:0.25 20:100:0.25 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.0p25.f000 | grep ':TMP:' |wgrib2 gdas.t00z.pgrb2.0p25.f000 -lola 100:100:0.25 20:100:0.25 out.bin bin -i
这样会生成一个含有所有温度信息的二进制文件。二进制文件的格式是
字节长度1 ,内容, 字节长度1,字节长度2, 内容, 字节长度2,……
内容直接就是float32类型的矩阵,由于我们导出的是一个区域,实际上有100x100个点。绘制脚本仅仅导出一部分数据进行绘制,否则太多了。
clear
lats = 100:0.25:124.75;
lons = 20:0.25:44.75;
[mx,my] = meshgrid(lons,lats);
fid = fopen('out.bin','rb');
asize = fread(fid,[1],'uint32');
allv = fread(fid,[100,100],'float32');
ct = 0;
while (asize==40000)
ct = ct + 1;
if (ct>=19 && ct <=27)
subplot(3,3,ct-18);
mesh(mx,my,allv);
endif
asize = fread(fid,[1],'uint32');
asize = fread(fid,[1],'uint32');
allv = fread(fid,[100,100],'float32');
endwhile
fclose(fid);
输出图像是几个海拔高度下的温度分布,单位为开尔文(K):
4 后记
在帮助学生完成了后续的可视化工作后,还是觉得意犹未尽。不知道学习气象的同学知不知道国内气象数据的开发情况,特别是有没有这样的开放数据网站。学生给我看的国内的天气预报接口都太花哨、却不够深入、方便——做地理信息的展示,原始数据是最合适的选择。
4.1关于数据网站
看到这个网站的美工,不由得想起90年代的门户。虽然很土,但是配合了如此的内容,真的有一种老派JazzBand的感觉。如果稍微深入阅读一下网站的说明,就会发现这个数据源背后,是一群从80年代就开始,用fortran写代码、在unix系统里做研究的超级神人,以及他们带出的学生。能够把如此多的科学参数,整合到一起,有序的管理起来,同时,开发出非常灵活、有效的工具链进行发布,本身就非常厉害了。更可贵的是,文档也非常通俗易懂,这种科学精神真的值得我们学习。
4.2留学生的计算机技能准备
出国留学的同学,一定会发现国外大学与国内的不同。有些大学论文要求用latex,周围同学做笔记,用的也是libreoffice等等开源文本处理器。计算实验用Matlab学生版是有的,也有VC++ 可用,但学习一下python, Octave,GNU C, Linux绝对没有坏处——这样想抄别人的也好抄。