grads绘图与编程(值得收藏).ppt

《grads绘图与编程(值得收藏).ppt》由会员分享，可在线阅读，更多相关《grads绘图与编程(值得收藏).ppt（39页珍藏版）》请在三一办公上搜索。

1、1,GrADS绘图与编程,2,图例1 1月份500hPa位势高度场,3,图例2 500hPa高度场（1、4、7、10月）,4,图例3 亚洲季风区7月850hPa流场图,5,图例4 SST和Sea Level分布图,6,图例5（130E，25N）处的感热和潜热的时间演变,7,第一章 概述,简介 GrADS的安装 GrADS软件的组成 GrADS的启动和退出 所需预备知识 学习本课程的要求,8,简 介,GrADS 是Grid Analysis and Display System的缩写,它由美国马里兰大学气象系Brian E.Doty开发。利用该系统可实现包括格点数据和站点数据的彩色气象图形分析和

2、显示。具有操作简便、功能丰富、图形美观、显示快速的特点。在国内外气象界得到广泛应用。版本：等,9,安 装,1、将系统软件拷入硬盘某一子目录下，如：c:grads 2、修改autoexec.bat文件，加上 path c:grads;%path%set gaddir=c:grads set gascrp=c:grads,10,GrADS 软件包的组成,Grads.exe 系统的核心文件,由此进入GrADS环境Dos4gw.exe 由GrADS.exe文件调用Gribmap.exe 产生格点资料映射文件(*.idx)Gribscan.exe 看Grib码资料，转为TXT格式文件Gxtran.exe

3、 显示*.gmf格式的图形文件Stnmap.exe 产生台站资料的映射文件(*.map)Wgrib.exe Grib码资料的解读程序Gv.exe 将*.gmf格式的图形文件转换为*.wmf 格式,11,启动和退出GrADS,在DOS下键入：grads既可进入；键入:quit即可退出grads,返回DOS。GrADS在初始化绘图环境之前，将提示用户选择风景画(Landscape)或肖像画(Portrait)形式。风景画形式的大小为118.5英寸，肖像画为 8.5 11英寸。启动GrADS系统时，有以下选择：-b 以批处理形式运行GrADS。-l 以风景画形式运行GrADS。-p 以肖像画形式运行

4、GrADS。-c 进入GrADS后首先执行随后跟在RUN命令后的描述文件，如：grads c“run profile.gs”上述选项可结合起来使用，如：grads blc“run batch.gs”,11,8.5,8.5,11,12,所需预备知识DOS和Windows基础 FORTRAN编程 WORD基础 学习要求强调上机编程实践,13,第二章 GrADS 数据格式,目标(1)能将文本格式数据转为GrADS格式(2)能写出数据描述文件(*.ctl)一、维数环境的概念,重要概念,1、含义：GrADS视每一个变量(VAR)场为一个四维数据集，即包括三维空间(x,y,z)和一维时间(t)。2、作用：

5、说明和指定随后的分析或图形操作时参加的原始数据集的维数范围。3、定义方法：Set lat|lon|lev|time val1 地球坐标 Set x|y|z|t val1 格点坐标,14,x(或lon)从西向东的水平坐标y(或lat)从南到北的水平坐标z(或lev)从地面到高空的垂直坐标t(或time)时序坐标。,如：set lon 40 160；set lat 0 60 定义了水平变化范围 set lat 30；set lon 0 180 定义了沿30N的纬向变化范围 set time jan81 dec94 定义了从81年1月到94年12月的时段所以，固定的维数环境和变化的维数环境相结合，就

6、构成了当前的维数环境。,15,(1)x,y,z,t均固定时，得到一个单值数据点(2)x,y,z,t中只有一维变化时，得到一维曲线(3)二维发生变化时，对应于二维剖面图，如X-Y平面图，X-Z，Y-Z，X-T，Y-T，Z-T剖面图(4)当三维发生变化时，GrADS以动画方式显示二维切片。,16,说出下列维数环境的含义,Set lon 60 150Set lat 0 30Set t 7Set lev 200,Set lon 120Set lat 30Set z 1 12Set t 1,Set lat 30Set lon 60 130Set t 1 12Set z 1,Set lon 60 150S

7、et lat 0 30Set t 1 12Set lev 200,17,二、文本格式数据 GrADS格式,现有如下资料：名称：U850，V850，U200，V200，H500 和TSFC范围：60150E，040N分辨率：2.52.5时间：1982年1月1985年12月的逐月资料。如何将上述资料写成GrADS下的数据格式？,18,顺序：(x,y)z v t,T=1T=2T=3,U VHT,U850(X,Y),U200(X,Y),时(v),分(z),秒(x,y),5月20日(t),Parameter(ii=37,jj=17)Real var(ii,jj)Open(1,file=u850.dat)

8、Open(2,file=v850.dat)Open(3,file=u200.dat)Open(4,file=v850.dat)Open(5,file=h500.dat)Open(6,file=tsfc.dat)Open(9,file=data.grd,form=unformatted,access=direct,recl=ii*jj*4)Irec=1Do 200 iy=1,4Do 100 m=1,12Read(1,1000)Read(1,2000)(var(I,j),I=1,ii),j=1,jj)Write(9,rec=irec)(var(I,j),I=1,ii),j=1,jj)Irec=ir

9、ec+1,文本格式数据源,GrADS格式数据,FORTRAN源程序,无格式，直接记录,Read(3,1000)Read(3,2000)(var(I,j),I=1,ii),j=1,jj)Write(9,rec=irec)(var(I,j)I=1,ii),j=1,jj)Irec=irec+1Read(2,1000)Read(2,2000)(var(I,j),I=1,ii),j=1,jj)Write(9,rec=irec)(var(I,j)I=1,ii),j=1,jj)Irec=irec+1Read(4,1000)Read(4,2000)(var(I,j),I=1,ii),j=1,jj)Write(

10、9,rec=irec)(var(I,j)I=1,ii),j=1,jj)Irec=irec+1Read(5,1000)Read(5,3000)(var(I,j),I=1,ii),j=1,jj)Write(9,rec=irec)(var(I,j)I=1,ii),j=1,jj)Irec=irec+1Read(6,1000)Read(6,4000)(var(I,j),I=1,ii),j=1,jj)Write(9,rec=irec)(var(I,j)I=1,ii),j=1,jj)Irec=irec+1,100 continue200 continue1000 format(2i7)2000 format

11、(37f6.2)3000 format(37f8.1)4000 format(37f7.2)end,22,三、数据描述文件,Data.ctl 文件内容：Dset data.grdUndef 9.99e+33Title Ncep/Ncar reanalysis projectXdef 37 linear 60 2.5Ydef 17 linear 0 2.5Zdef 2 levels 850 200Tdef 48 linear jan1982 1moVars 4U 2 99 u wind(m/s)V 2 99 v wind(m/s)H 1 99 H500T 1 99 Tsfc dataendvar

12、s,关键词：如dset,undef,title,xdef，ydef 等被描述的数据文件名缺测标记标题,确定维数环境,指定变量名,起始值间 隔层次数,23,数据描述文件一般包含以下几项：（1）被描述的数据文件名（DSET）（2）该数据说明文件的标题（TITLE）（3）数据类型、格式和选项（DTYPE，FORMAT，OPTION）（4）时间、空间维数环境设置（XDEF，YDEF，ZDEF，TDEF）（5）变量定义（VARS，ENDVARS）以下详细说明数据描述文件中各记录的含义：1、DSET data-set-name 给定二进制原始数据文件的文件名（可包含路径），若数据文件与描述文件在同一路径下

13、，可用省缺路径符号“”代表。,24,2、title string用字符串string简略描述数据文件的内容。3、undef value定义缺测值，GrADS在运算操作时和图形操作时将忽略这些格点。4、options 可以是其中：sequential表示数据是以顺序无格式形式存放，每个记录为一个(x,y)场;yrev表示 y维数方向与 ydef中说明的方向相反，即为由北到南；zrev表示 z维数方向与 zdef中说明的方向相反，即为由 上到下；byteswapped表示二进制数据的位存放顺序取反序；big_endian和 little_endian不同机器之间二进制位存放顺序的自动改变；temp

14、late表示 多个时间序列原始数据文件共用一个数据描述文件。例如：,25,有1may92.dat,2may92.dat,.等数据文件，需用一个共同的描述文件，这时的CTL文件可写为：Dset%d1%mc%y2.dat.Options templateTdef 72 linear 0z1may92 1 dy其中：%y2:2位数年；%y4:4位数年；%m1:1或2位数月；%m2:2位数月；%mc:3个字符的月份英文缩写;%d1：1或2位数天；%d2：2位数天；%h1：1或2位数小时；%h2：2位数小时5、xdef number linear start inc 或 xdef number leve

15、ls 设置网格点与经度的对应关系。其中number是x方向网格点数，linear或 levels表明网格点映射类型，start起始经度，inc格距大小，value-list表示X方向各格点的列表。,26,6、ydef number mapping start inc 或 ydef number levels values-list设置Y方向格点与纬度的映射关系。其中number为Y方向的格点数，mapping表示映射方式，有：Linear:线性映射Gausr15：高斯R15纬度，Gausr20,gausr30,gausr40等对于线性映射linear，start 为起始纬度，inc为Y方向格距

16、。对于高斯映射，start为第一高斯网格数。对于levels映射，value-list为Y方向取值表7、zdef number linear start inc 或 zdef number levels value-list设置气压面与垂直网格点的映射关系。如：Zdef 10 linear 1000 100Zdef 10 levels 1000,925,850,700,600,500,400,300,200,100,27,8、tdef number linear start inc设置网格值与时间的映射关系。其中：number表示数据集中的时次数，start表示起始时间/日期，用GrADS绝对

17、时间表示，其格式为：Hh:mm Z dd mmm yyyyInc表示时间增量，其格式为VVKK，（VV为增量值，KK增量类型，有mm分，hr小时，dy天，mo月，yr年)9、vars number（变量个数）变量名 垂直层次数 预留值 变量说明 endvars,28,地形高度资料数据描述文件,Dset orog.grdUndef 9.99e+33Title Ncep/Ncar reanalysis projectXdef 37 linear 60 2.5Ydef 17 linear 0 2.5Zdef 1 levels 850Tdef 1 linear jan1982 1 moVars 1Or

18、og 0 99 orographyendvars,29,第三章 GrADS绘图和参数设置,30,以BAR形式画出沿35N的地形高度剖面图（基点取为3000米）,open orog.ctlc;set grads offset lat 35set xlopts 1 5 0.16;set ylopts 1 5 0.16set gxout bar;set bargap 50set barbase 3000d orogq w2xy 60 3000X1=subwrd(result,3);y1=subwrd(result,6)q w2xy 150 3000X2=subwrd(result,3);y2=sub

19、wrd(result,6)set line 2 1 6draw line x1 y1 x2 y2 draw title OROG ALONG 35Ndraw xlab lon;draw ylab OROG(M),画水平轴线,31,画出地形高度等值线图（间隔500M，等值线标注为500M字样，大于3000M的区域打阴影）,open orogset parea 1 10 1 8c;set grid off;set grads offset xlopts 1 5 0.16;set ylopts 1 5 0.16set map 1 1 10set gxout shaded;set rbcols 5se

20、t cmin 3000;d orogset gxout contour;set cint 500set clab%gm;d orogdraw title OROGdraw xlab lon;draw ylab lat,32,画出1、7月份的H500气候场,open dataset parea 1 10 1 8c;set grads off;set grid offset xlopts 1 5 0.16;set ylopts 1 5 0.16set map 1 1 10define h1pj=ave(h,t=1,t=48,12)define h7pj=ave(h,t=7,t=48,12)d h1

21、pj;draw title H500 in Jandraw xlab lon;draw ylab latc;set grads off;d h7pjdraw title H500 in Juldraw xlab lon;draw ylab lat,33,画出84年1月、7月的H500距平场（负值区加阴影）,open dataset parea 1 10 1 8;set grads off;set grid offset xlopts 1 5 0.16;set ylopts 1 5 0.16;set map 1 1 10define h1pj=ave(h,t=1,t=48,12);define

22、h7pj=ave(h,t=7,t=48,12)define h84jan=h(time=jan84)-h1pj;define h84jul=h(time=jul84)-h7pjset gxout shaded;set rbcols 5;set cmax o;d h84janset gxout contour;d h84jandraw clevs 0;set cthick 10;d h84jandraw title H500 anomaly in jan 1984;draw xlab lon;set ylab latc;set grads off;set gxout shaded;set rbc

23、ols 5;set cmax 0;d h84julset gxout contour;set cint 5;d h84julset clevs 0;set cthick 10;d h84juldraw title H500 anomaly in Jul 1984;draw xlab lon;set ylab lat,34,画出UV850的1月份气候流场图（在气旋和反气旋中心位置标注字符C和A，其中地形高度大于1500M的区域画阴影）,open data;open orogset parea 1 10 1 8;set grads off;set grid offset xlopts 1 5 0.

24、16;set ylopts 1 5 0.16;set map 1 1 10define upj=ave(u,t=1,t=48,12);define vpj=ave(v,t=1,t=48,12)define uupj=maskout(upj,1500-orog.2)define vvpj=maskout(vpj,1500-orog.2)set gxout stream;set strmden 4;d uupj;vvpjset gxout contour;set cthick 10;set clevs 1500;d orog.2q w2xy 88 35;x2=subwrd(result,3);y2

25、=subwrd(result,6)set string 2 c 6 0;set strsiz 0.2;draw string x2 y2 TIBETAN PLATEAUq w2xy 82 20;x1=subwrd(result,3);y1=subwrd(result,6)set string 2 c 6 0;set strsiz 0.2;draw string x1 y1 Aq w2xy 119 25;x3=subwrd(result,3);y3=subwrd(result,6)set string 2 c 6 0;set strsiz 0.2;draw string x3 y3 Adraw

26、title UV850 flowfield in Jandraw xlab lon;draw ylab lat,35,画出82年3月、83年4月、85年7月合成的H500场,open dataset parea 1 10 1 8set xlopts 1 5 0.16;set ylopts 1 5 0.16set map 1 1 10c;set grads off;set grid offdefine aa=(h(time=mar82)+h(time=apr83)+h(time=jul85)/3d aadraw title H500 for 82.3/83.4/85.7draw xlab lon

27、;draw ylab lat,36,画出H500单点相关图相关基点为（100120E，2030N）平均,open data;set t 1 48define one=aave(h,lon=110,lon=120,lat=20,lat=30)set t 1define hbar=ave(h,t=1,t=48)define onebar=ave(one,t=1,t=48)define hdd=sqrt(ave(pow(hbar-h),2),t=1,t=48)define onedd=sqrt(ave(pow(onebar-one,2),t=1,t=48)define hone=ave(h-hbar

28、)*(one-onebar),t=1,t=48)define r=hone/(hdd*onedd)d smth9(r)draw title H500 one point corr with base point at 110-120E,20-30N,37,逐月显示82年1月85年12月的UV850风矢量图,open data;set map 2 1 10;set xlopts 1 5 0.16;set ylopts 1 5 0.16set grid offIy=82While(iy=85)Mon=1While(mon=12)Tt=(iy-82)*12+monset t tt d u;vdraw

29、 title UV850 for iy.monPull dummy;cMon=mon+1EndwhileIy=iy+1endwhile,二重循环,38,画出（60100E，510N）平均的U850-U200的逐月变化曲线在大于30和小于-15的峰谷值点位置画,open data;set lon 120;set lat 20;set t 1 48define wb=aave(u(lev=850)-u(lev=200),lon=60,lon=100,lat=5,lat=20)d wbP=1While(P30|r-15)q gr2xy P r X=subwrd(result,3);y=subwrd(

30、result,6)Draw mark 3 x y 0.15EndifP=P+1endwhile,39,在北半球极地投影下标经度值,*Draw lon values for npj projection figure*to run this function,just in your*.ps add:*lon 0 30 60 90 120 150 180 210 240 270 300 330Function abc(arg)I=1;while(I0);draw string x1 y1 L.IE;endifIf(L.I=180);draw string x1 y1 180;endifIf(L.I180);lll=360-L.I;draw string x1 y1 lllW;endifI=I+1Endwhilereturn,文件名LON.GS,假设最多标12个经度值，读入欲标经度值，放入复合变量中,标值位置离最外纬圈3度,