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

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开发。利用该系统可实现包括格点数据和站点数据的彩色气象图形分析和显示。具有操作简便、功能丰富、图形美观、显示快速的特点。在国内外气象界得到广泛应用。版本：等,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 显示*.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 以肖像画形式运行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、作用：说明和指定随后的分析或图形操作时参加的原始数据集的维数范围。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月的时段所以，固定的维数环境和变化的维数环境相结合，就构成了当前的维数环境。,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 150Set 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)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=irec+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(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(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 dataendvars,关键词：如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 给定二进制原始数据文件的文件名（可包含路径），若数据文件与描述文件在同一路径下，可用省缺路径符号“”代表。,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不同机器之间二进制位存放顺序的自动改变；template表示 多个时间序列原始数据文件共用一个数据描述文件。例如：,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 levels 设置网格点与经度的对应关系。其中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方向格距。对于高斯映射，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绝对时间表示，其格式为：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 1Orog 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=subwrd(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 5set 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 h1pj;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 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 rbcols 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.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=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 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;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)*(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 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(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度,