很多出國深造的同學，都對國外高校中的計算機教學、使用記憶猶新。國内一般院校的老師很多都是從微軟的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

httpd服務的部署、啟動、配置和簡單優化一、部署二、啟動三、配置檔案

配置網頁内容通路

手動安裝Intel network I217-LM網卡的Linux驅動

禁止ubuntu系統彈出報錯界面

Ubuntu Linux下Apache的配置檔案

samba伺服器的功能

【Linux】UDP廣播封包接收速率問題

Linux裝置模型（中）之上層容器

PowerPC平台 Linux移植三