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

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

2.1各類數(shù)據(jù)格式的特點一般二進制無格式直接或者順序記錄的數(shù)據(jù)格式：可以是格點數(shù)據(jù)或站點數(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ù)文件的轉換

1.轉換方法

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

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

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

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

y（緯度）：從南到北

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

vars（各種物理變量）

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

任何一維可省略。

12345………..x

y………54321Z3.舉例（格點資料）

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

編程時數(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時次：月資料，4年，共48個月Fortran程序編寫如下：

!定義一個X,Y,Z方向的格點數(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程序在一個文件夾下，打開數(shù)據(jù)文件時可以不寫路徑。

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

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重新進行記錄累加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')!打開目標文件open(12,file='mhy.grd',form='binary')

!打開原始數(shù)據(jù)文件,如果數(shù)據(jù)文件和Fortran程序在一個文件!夾下，打開數(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ù)文件讀入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ù)資料寫入新文件中注意事項：二進制寫入的數(shù)據(jù)必須是實型的(real)若數(shù)據(jù)為整型，寫入時要轉換成實型的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ù)結構、變量描述等等。2.數(shù)據(jù)描述文件的形式

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

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

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

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

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

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

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

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

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

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

LINEAR指明是等間隔分布格點，

Start起點坐標，負數(shù)表示西經(jīng)；

increment網(wǎng)格間距。

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

設置Y方向格點與緯度的對應關系YDEFnumber

<LINEARStartincrement>或YDEFnumber

<LEVELSvalue-list>或YDEFnumber<mapping

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

表示共有20個Y方向網(wǎng)格點，起始點為高斯R40網(wǎng)格下的高斯緯度15（即64.100S）,實際這20個網(wǎ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設置垂直網(wǎng)格點與氣壓面的對應關系ZDEFnumber<LINEARStartincrement>或ZDEFnumber<LEVELSvalue-list>

zdef7levels1000850700500300200100

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

TDEFnumberLINEARstart-timeincrement

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

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

分鐘6hr

小時1dy

天1mo

月1yr

年例如：12Z1JAN199014:20Z22JAN1987JUN1960

Z：世界時12Z代表北京時間20點VARSnumber

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

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

dtype用于說明數(shù)據(jù)格式的類型當數(shù)據(jù)為格點binary格式時省略dtype的可選參數(shù)

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

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

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

注意事項：如果電腦中常用文件類型隱藏時，請修改該設置。打開“我的電腦”，打開在窗口第一行第五個“工具”，選“文件夾選項”，點“查看”，將滾動條下拉至“隱藏已知文件類型的擴展名”，將其前面對號去掉，點“確定”即可。鼠標右鍵新建文本文檔，打開，寫入ctl內容，然后保存為file.ctl。練習：為mhy.grd文件（剛才進行數(shù)據(jù)轉換得到的新文件）編寫數(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ù)內容說明(自描述)一般分辨率為2.5°×2.5°或者1°×1°

標準NC格式,GrADS可直接處理，fortranga->sdfopenfilename.nc?特別提醒：雖然GrADS支持這種格式的數(shù)據(jù)，但是涉及到要用程序對數(shù)據(jù)進行處理時，F(xiàn)ortran程序并不支持這種格式的數(shù)據(jù)。如果只是使用這類數(shù)據(jù)進行畫圖和GrADS的命令處理，則可以直接操作。如果要用于程序使用的話，應該利用命令“fwrite”將所用數(shù)據(jù)提取為“.grd”的二進制數(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、從“開始”菜單的“運行”選項里輸入“cmd”進入DOS狀態(tài)，然后cdc:\PCGrADS\win32e。2、在該路徑下輸入命令：

grib2ctl路徑/model.grb

>

路徑/model.ctl

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

gribmap-i路徑/model.ctl

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

open路徑/model.ctl

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

年份月份日期臺站號緯度經(jīng)度海拔高度日平均氣壓(0.1hpa)日平均氣溫(0.1℃)……..1951-2005年7月160站月降水量160站站點號和經(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)！讀入站點號和經(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對數(shù)據(jù)進行處理！距平值寫成grads可以識別的站點數(shù)據(jù)！將整數(shù)數(shù)組idx轉化成字符型數(shù)組id，id是8個字符，這個是固定的

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

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

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

ilev=0！取0表示一個時次記錄的結束

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

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