第2章 數(shù)據(jù)處理

第2章數(shù)據(jù)處理在使用GrADS繪圖前，必須使所用數(shù)據(jù)文件滿足GrADS的數(shù)據(jù)格式要求。

2.1各類數(shù)據(jù)格式的特點(diǎn)一般二進(jìn)制無格式直接或者順序記錄的數(shù)據(jù)格式：可以是格點(diǎn)數(shù)據(jù)或站點(diǎn)數(shù)據(jù)，可以用Fortran或者C語言讀寫。需要數(shù)據(jù)描述文件.ctl

GRIB碼數(shù)據(jù)格式：是自定義的數(shù)據(jù)格式，壓縮率高，占用空間小。例如NCEP提供的數(shù)據(jù)資料。NETCDF（.nc）等通用數(shù)據(jù)格式：自定義數(shù)據(jù)格式，精確性好，便于傳輸。

文件中自帶描述文件。2.2數(shù)據(jù)文件的轉(zhuǎn)換

1.轉(zhuǎn)換方法

轉(zhuǎn)換文件的數(shù)據(jù)存放格式，一般是通過PowerStation或VisualFortran等軟件使用Fortran或者C語言來編程轉(zhuǎn)換。

2.數(shù)據(jù)存放形式

二進(jìn)制數(shù)據(jù)排放順序從內(nèi)循環(huán)到外循環(huán)依次是：

x（經(jīng)度）：從西到東

y（緯度）：從南到北

z（高度層數(shù)）：從低層到高層

vars（各種物理變量）

t（時(shí)次）x(lon)y(lat)z(lev)vars(不同變量)time

任何一維可省略。

12345………..x

y………54321Z3.舉例（格點(diǎn)資料）

現(xiàn)有ASCII碼（十進(jìn)制存放格式）數(shù)據(jù)資料文件u.dat、v.dat和sst.dat，其空間范圍：60-150°E，0-40°N；層次：u、v為850、200hPa；時(shí)段：1982.1-1985.12；分辨率：2.5*2.5。要求編寫出將這三個(gè)文件轉(zhuǎn)換成二進(jìn)制無格式直接存?。℅rads格式.grd）文件的Fortran程序。

編程時(shí)數(shù)組大?。篨方向：150-60=90，90/2.5+1=37Y方向：40-0=40，40/2.5+1=17Z方向：u、v為850、200hPa

，Z為2sst為海表溫度，只有一層，Z為1T時(shí)次：月資料，4年，共48個(gè)月Fortran程序編寫如下：

!定義一個(gè)X,Y,Z方向的格點(diǎn)數(shù)以及總時(shí)次ntparameter(nx=37,ny=17,nz=2,nt=48)!定義數(shù)組

dimensionu(nx,ny,nz,nt),v(nx,ny,nz,nt),sst(nx,ny,nt)

!打開原始數(shù)據(jù)文件,如果數(shù)據(jù)文件和Fortran程序在一個(gè)文件夾下，打開數(shù)據(jù)文件時(shí)可以不寫路徑。

open(1,file='u.dat')open(2,file='v.dat')open(3,file='sst.dat')!打開目標(biāo)文件，recl為記錄的長度，對(duì)于直接存取文件必須指定記錄長度。

open(12,file='mhy.grd',form='unformatted',#access='direct',recl=nx*ny*4)

!把數(shù)據(jù)文件讀入do100it=1,ntdo50iz=1,nzread(1,*)((u(i,j,iz,it),i=1,nx),j=1,ny)50continuedo51iz=1,nzread(2,*)((v(i,j,iz,it),i=1,nx),j=1,ny)51continueread(3,*)((sst(i,j,it),i=1,nx),j=1,ny)100continue!將數(shù)據(jù)資料寫入新文件中，irec重新進(jìn)行記錄累加irec=0do200it=1,ntdo30iz=1,nzirec=irec+1write(12,rec=irec)((u(i,j,iz,it),i=1,nx),j=1,ny)30Continuedo31iz=1,nzirec=irec+1write(12,rec=irec)((v(i,j,iz,it),i=1,nx),j=1,ny)31continueirec=irec+1write(12,rec=irec)((sst(i,j,it),i=1,nx),j=1,ny)200continueendirec=0do200it=1,ntdo30iz=1,nzirec=irec+1write(12,rec=irec)((u(i,j,iz,it),i=1,nx),j=1,ny)30Continuedo31iz=1,nzirec=irec+1write(12,rec=irec)((v(i,j,iz,it),i=1,nx),j=1,ny)31continueirec=irec+1write(12,rec=irec)((sst(i,j,it),i=1,nx),j=1,ny)200continue如果不考慮采用無格式直接方式存取文件的話，可以使用下面的編程形式，較為簡單。parameter(nx=37,ny=17,nz=2,nt=48)dimensionu(nx,ny,nz,nt),v(nx,ny,nz,nt),sst(nx,ny,nt)open(1,file='u.dat')open(2,file='v.dat')open(3,file='sst.dat')!打開目標(biāo)文件open(12,file='mhy.grd',form='binary')

!打開原始數(shù)據(jù)文件,如果數(shù)據(jù)文件和Fortran程序在一個(gè)文件!夾下，打開數(shù)據(jù)文件時(shí)可以不寫路徑do100it=1,ntdo50iz=1,nzread(1,*)((u(i,j,iz,it),i=1,nx),j=1,ny)50continuedo51iz=1,nzread(2,*)((v(i,j,iz,it),i=1,nx),j=1,ny)51continueread(3,*)((sst(i,j,it),i=1,nx),j=1,ny)100continue!把數(shù)據(jù)文件讀入do200it=1,ntdo30iz=1,nzwrite(12)((u(i,j,iz,it),i=1,nx),j=1,ny)30Continuedo31iz=1,nzwrite(12)((v(i,j,iz,it),i=1,nx),j=1,ny)31continuewrite(12)((sst(i,j,it),i=1,nx),j=1,ny)200continueend!將數(shù)據(jù)資料寫入新文件中注意事項(xiàng)：二進(jìn)制寫入的數(shù)據(jù)必須是實(shí)型的(real)若數(shù)據(jù)為整型，寫入時(shí)要轉(zhuǎn)換成實(shí)型的real(iu(i,j,k,it))2.2數(shù)據(jù)描述文件（.ctl）

GrADS繪圖不能直接使用“數(shù)據(jù)文件”，而是通過“數(shù)據(jù)描述文件”間接使用“數(shù)據(jù)文件”。1.數(shù)據(jù)描述文件的概念

數(shù)據(jù)描述文件是原始數(shù)據(jù)文件的描述文件。用以描述原始數(shù)據(jù)集的基本信息，包括數(shù)據(jù)集文件名、數(shù)據(jù)類型、數(shù)據(jù)結(jié)構(gòu)、變量描述等等。2.數(shù)據(jù)描述文件的形式

OPTIONS<keywords>定義數(shù)據(jù)格式選項(xiàng)，keywords有：<yrev><zrev><sequential><byteswapped><template><big-endian><little-endian>

缺省時(shí)為direct直接記錄格式（form='binary'）<yrev>自北向南的數(shù)據(jù)<zrev>自高層到低層的數(shù)據(jù)<sequential>順序無格式存貯的數(shù)據(jù)<byteswapped>二進(jìn)制數(shù)據(jù)的位存放順序取反序

<big-endian>在SUN,SGI等機(jī)器上生成的數(shù)據(jù)<little-endian>在PC,DEC等機(jī)器上生成的數(shù)據(jù)服務(wù)器個(gè)人電腦big-endian、little-endian用于自動(dòng)改變二進(jìn)制位存放順序<template>

多個(gè)時(shí)間序列原始數(shù)據(jù)文件用一個(gè)數(shù)據(jù)描述文件統(tǒng)一地描述這些原始數(shù)據(jù)時(shí)采用的選項(xiàng)，這些數(shù)據(jù)文件的文件名形式由dset定義的形式命名文件名，提示所含數(shù)據(jù)的時(shí)次。例如：一個(gè)逐日的數(shù)據(jù)集每天數(shù)據(jù)放到一個(gè)文件中，每個(gè)文件名形式為：

1may92.dat，2may92.dat，3may92.datdset%d1%mc%y2.datoptionstemplatetdef30linear0Z1may19921dy

定義時(shí)間范圍和增量例如：正確的替換為：%y2兩位數(shù)年%y44位數(shù)年

%m11或2位數(shù)月%m22位數(shù)月（用0補(bǔ)齊1位數(shù)）

%mc3字符月份縮寫%d11或2位數(shù)天%d22位數(shù)天

%h11或2位數(shù)小時(shí)%h22位小時(shí)

設(shè)置X方向格點(diǎn)與經(jīng)度的對(duì)應(yīng)關(guān)系XDEFnumber<LINEARStartincrement>或XDEFnumber<LEVELSvalue-list>number（>=1）給定格點(diǎn)數(shù)，整形；

LINEAR指明是等間隔分布格點(diǎn)，

Start起點(diǎn)坐標(biāo)，負(fù)數(shù)表示西經(jīng)；

increment網(wǎng)格間距。

LEVELS參數(shù)指明是不等間隔分布格點(diǎn)，列出具體每個(gè)格點(diǎn)的坐標(biāo)值（以空格分開）。

設(shè)置Y方向格點(diǎn)與緯度的對(duì)應(yīng)關(guān)系YDEFnumber

<LINEARStartincrement>或YDEFnumber

<LEVELSvalue-list>或YDEFnumber<mapping

Start>mapping表示映射方式，有以下幾種：GAUSR15--高斯R15緯度(Gaussian)GAUSR20--高斯R20緯度GAUSR30--高斯R30緯度GAUSR40--高斯R40緯度在linear中，Start起點(diǎn)坐標(biāo)負(fù)數(shù)表示南緯；在高斯映射中，Start起點(diǎn)坐標(biāo)表示高斯緯度。方向:自西向東,自南向北(永遠(yuǎn)不變！)單位:個(gè)緯度/經(jīng)度正負(fù)值:以東為正/以北為正xdef144linear02.5ydef73linear-902.5方向維數(shù)定義格點(diǎn)數(shù)坐標(biāo)值線性增加起始坐標(biāo)坐標(biāo)間隔YDEF20GAUSR4015

表示共有20個(gè)Y方向網(wǎng)格點(diǎn)，起始點(diǎn)為高斯R40網(wǎng)格下的高斯緯度15（即64.100S）,實(shí)際這20個(gè)網(wǎng)格點(diǎn)對(duì)應(yīng)的緯度值為：-64.10,-62.34,-60.58,-58.83,-57.07,-55.32,-53.56,-51.80,-50.05,-48.29,-46.54,-44.78,-43.02,-41.27,-39.51,-37.76,-36.00,-34.24,-32.49,-30.73設(shè)置垂直網(wǎng)格點(diǎn)與氣壓面的對(duì)應(yīng)關(guān)系ZDEFnumber<LINEARStartincrement>或ZDEFnumber<LEVELSvalue-list>

zdef7levels1000850700500300200100

列出具體坐標(biāo)值方向:自下而上單位:層hPa/其他

TDEFnumberLINEARstart-timeincrement

設(shè)置網(wǎng)格值與時(shí)間的映射關(guān)系。時(shí)間表示法：hh:mmZddmmmyyyyhourminutedaymonyeartdef12linear00:00z21jan20001mo

時(shí)間步長只有l(wèi)inear關(guān)鍵字起始時(shí)間必須按照上述格式寫年份用4位數(shù)寫21日00:00時(shí)按情況可省略2000年1月(不可省略)3mn

分鐘6hr

小時(shí)1dy

天1mo

月1yr

年例如：12Z1JAN199014:20Z22JAN1987JUN1960

Z：世界時(shí)12Z代表北京時(shí)間20點(diǎn)VARSnumber

變量描述開始，同時(shí)給出變量個(gè)數(shù)number。

abrevlevsunitsdescription………..abrev，變量名稱，為1到8個(gè)字符組成的該變量的縮寫名，畫圖或處理數(shù)據(jù)將用到，以字母開頭levs，數(shù)字，代表變量層數(shù)，0表示只有一層units，單位，為GRIB預(yù)留，給99description，對(duì)變量的文字描述，最多40個(gè)字符用ENDVARS表示數(shù)據(jù)描述文件結(jié)束。

dtype用于說明數(shù)據(jù)格式的類型當(dāng)數(shù)據(jù)為格點(diǎn)binary格式時(shí)省略dtype的可選參數(shù)

–grib–station站點(diǎn)數(shù)據(jù)

–hdfsds–netcdfv1.9版本新增?ctl書寫注意事項(xiàng)???每行頂格寫;雖然大小寫沒關(guān)系,推薦在一個(gè)關(guān)鍵字中都用大寫或都用小寫;如寫成dset,而不寫成DsEt；每行記錄的各項(xiàng)以空格分開；每行記錄不超過80個(gè)字符；注釋行請(qǐng)用*開始,并頂格寫;不能出現(xiàn)在變量列表中，即：在vars~endvars之間不要寫注釋行;不要在GrADS中使用制表(tab)鍵!4.數(shù)據(jù)描述文件的建立

數(shù)據(jù)描述文件為純文本格式，可用一般的編輯器產(chǎn)生（如記事本，寫字板等）。描述文件編輯完成后，保存時(shí)為.ctl。

注意事項(xiàng)：如果電腦中常用文件類型隱藏時(shí)，請(qǐng)修改該設(shè)置。打開“我的電腦”，打開在窗口第一行第五個(gè)“工具”，選“文件夾選項(xiàng)”，點(diǎn)“查看”，將滾動(dòng)條下拉至“隱藏已知文件類型的擴(kuò)展名”，將其前面對(duì)號(hào)去掉，點(diǎn)“確定”即可。鼠標(biāo)右鍵新建文本文檔，打開，寫入ctl內(nèi)容，然后保存為file.ctl。練習(xí)：為mhy.grd文件（剛才進(jìn)行數(shù)據(jù)轉(zhuǎn)換得到的新文件）編寫數(shù)據(jù)描述文件mhy.ctl。dsetmhy.grdundef-9.99E+33titleNCEP/NCARREANALYSISPROJECTxdef37linear60.02.5ydef17linear0.02.5zdef2levels850200tdef48linearJAN19821movars3u299uwind(m/s)v299vwind(m/s)sst099TSFCDATAendvars打開ctlopenh:\grads\sst.win.ctl5、netCDF格式數(shù)據(jù)–ECMWF/NCEP再分析資料等–常用后綴.nc–自帶數(shù)據(jù)內(nèi)容說明(自描述)一般分辨率為2.5°×2.5°或者1°×1°

標(biāo)準(zhǔn)NC格式,GrADS可直接處理，fortranga->sdfopenfilename.nc?特別提醒：雖然GrADS支持這種格式的數(shù)據(jù)，但是涉及到要用程序?qū)?shù)據(jù)進(jìn)行處理時(shí)，F(xiàn)ortran程序并不支持這種格式的數(shù)據(jù)。如果只是使用這類數(shù)據(jù)進(jìn)行畫圖和GrADS的命令處理，則可以直接操作。如果要用于程序使用的話，應(yīng)該利用命令“fwrite”將所用數(shù)據(jù)提取為“.grd”的二進(jìn)制數(shù)據(jù)文件，以后再用于編程。6、GRIB格式數(shù)據(jù)T213,NECP/NCAR等常用后綴.grb需.ctl

和.idx

文件一起使用用grib2ctl.exe生成數(shù)據(jù)說明(描述)文件.ctl用gribmap–itest.ctl生成映射文件.idx1.8版本需另外下載放到win32e下！用grib2ctl.pl生成.ctl的方法：1、從“開始”菜單的“運(yùn)行”選項(xiàng)里輸入“cmd”進(jìn)入DOS狀態(tài)，然后cdc:\PCGrADS\win32e。2、在該路徑下輸入命令：

grib2ctl路徑/model.grb

>

路徑/model.ctl

這樣就可以在指定路徑下生成數(shù)據(jù)文件model.grb的數(shù)據(jù)描述文件model.ctl.注意：命令中間的“>”符號(hào)不能缺少，生成的數(shù)據(jù)描述文件的路徑和文件名都是用戶自己指定的。Grib2ctl.exe所在文件夾下用gribmap–itest.ctl生成映射文件.idx：繼續(xù)在C:\Program\PCGrADS\win32e下輸入命令：

gribmap-i路徑/model.ctl

這樣在原來存放數(shù)據(jù)文件model.grd的目錄下自動(dòng)生成數(shù)據(jù)文件的映射文件“model.idx”。在GrADS的命令窗口ga->提示符后輸入命令：

open路徑/model.ctl

即可使用model.grd數(shù)據(jù).GRIB格式數(shù)據(jù)ctl寫法特別提醒：雖然GrADS支持grib格式的數(shù)據(jù)，但是涉及到要用程序?qū)?shù)據(jù)進(jìn)行處理時(shí)，F(xiàn)ortran程序并不支持這種格式的數(shù)據(jù)。如果只是使用這類數(shù)據(jù)進(jìn)行畫圖和GrADS的命令處理，則可以直接操作。如果要用于程序使用的話，應(yīng)該利用命令“fwrite”將所用數(shù)據(jù)提取為“.grd”的二進(jìn)制數(shù)據(jù)文件，以后再用于編程。7、站點(diǎn)數(shù)據(jù)文本格式，點(diǎn)開可見數(shù)據(jù)值50136臺(tái)站06-07年逐日要素值，有說明文件

年份月份日期臺(tái)站號(hào)緯度經(jīng)度海拔高度日平均氣壓(0.1hpa)日平均氣溫(0.1℃)……..1951-2005年7月160站月降水量160站站點(diǎn)號(hào)和經(jīng)緯度文件：

idlonlat50353 51.72 126.6550632 48.77 121.9250527 49.22 119.7550434 50.45 121.7050557 49.17 125.2350745 47.38 123.9250756 47.43 126.9750788 47.23 131.9850873 46.49 130.1750978 45.28 130.9550953 45.75 126.7754094 44.57 129.60 dimensionrain(160,55) dimensionhp(160),ae(160,55) dimensionalat(160),alon(160)integeridx(160) character*8id(160)！讀入站點(diǎn)號(hào)和經(jīng)緯度資料

open(21,file='id160.dat')doist=1,160read(21,'(i5,2f7.2)')idx(ist),alat(ist),alon(ist)enddoclose(21)！讀入降水資料

open(41,file='r1607')doiy=1,55read(41,*)(rain(ii,iy),ii=1,160)enddoclose(41)!求55年平均值

doidx=1,160hp(idx)=0.0doit=1,55hp(idx)=hp(idx)+rain(idx,it)enddohp(idx)=hp(idx)/55enddo!求距平值

doidx=1,160doit=1,55ae(idx,it)=rain(idx,it)-hp(idx)enddoenddo對(duì)數(shù)據(jù)進(jìn)行處理！距平值寫成grads可以識(shí)別的站點(diǎn)數(shù)據(jù)！將整數(shù)數(shù)組idx轉(zhuǎn)化成字符型數(shù)組id，id是8個(gè)字符，這個(gè)是固定的

doi=1,160write(id,‘(i5)’)idx(i)

!或者id(i)=char(idx(i))enddo！打開新文件，存儲(chǔ)站點(diǎn)數(shù)據(jù)open(41,file=‘a(chǎn)rain7.tu.grd',form='binary')doit=1,55iflag=1！iflag”取1表示有地面變量，取0表示記錄中沒有地面觀測(cè)變量。ilev=1！ilev取1表示一個(gè)時(shí)次記錄的開始，取0表示一個(gè)時(shí)次記錄的結(jié)束time=0.0！time為本報(bào)的時(shí)次，不是一個(gè)確定的時(shí)刻值，只是該時(shí)次的！相對(duì)值

doist=1,160write(41)id(ist),alat(ist),alon(ist),time,ilev,iflag,ae(ist,it)enddo！給出一個(gè)特別記錄，表示這個(gè)時(shí)間組的記錄結(jié)束

ilev=0！取0表示一個(gè)時(shí)次記錄的結(jié)束

write(41)id(160),alat(160),alon(160),time,ilev,iflagenddoclose(41)enddsetd:\data\arain7.tu.grddtypestationstnmapd:\data\rain.mapundef-999.0titleraintdef55linearJUl19511yrvars1rain099anomalyrainfalldataENDVARS站點(diǎn)數(shù)據(jù)Ctl:站點(diǎn)映射文件名字stnmap用于將觀測(cè)站點(diǎn)數(shù)據(jù)與圖形顯示時(shí)的站點(diǎn)圖建立鏈接以加速顯示。語法為：

ga->stnmap-i站點(diǎn)數(shù)據(jù)描述文件路徑/**.ctl也可直接進(jìn)入stnmap交互式操作：stnmap.exe位于pcgrads目錄中win