《GrADS基础绘》PPT课件.ppt

GrADS基础绘图,何晏春2008-2-23,提纲,1、grads功能概述（1）、软件性能介绍:主要是与matlab对比（2）、grads在win/Linux下的下载、安装、运行：linux/windows环境变量的设置，安装后各个目录文件的内容（3）、紧接着进行grads的简单操作：对nc资料的操作，画两张简单的图：contour与shaded叠加（2）、Grads网络资源 2、Grads的图形显示类型及绘图设置（1）、一维：直方、折现、散点、误差图.（2）、二维：等值线、阴影图、风矢量、流线图、站点天气图（3）、图形设置：坐标、线型、页面、地图投影、地表状况、图形信息开关 3,Grads 函数 4、grads 批处理脚本 3、Grads数据准备（1）、ncep资料的直接读取（2）、二进制格点数据的生成、数据描述文件ctl作用 介绍grads格点数据的存储结构，即维数环境，生成二进制的fortran程序段；ctl文件的作用，ctl文件的结构。（3）、mm5out_to_Grads,wrfoutput_to_grads,（4）、站点数据的准备以及站点图形的绘制 4、grads函数的使用：（1）、数学函数(2）、统计函数（3）、网格处理函数.5、实例示范（1）precluster应用（2）pc演示，用grads的totourial。,参考资料,1、GrADS Tutorial2、GrADS Reference Card3、LASG:GrADS Mannual4、LASG:GrADS 实用绘图讲义，funny动力论坛5、动力论坛GrADS中文讨论区,常见软件,Panoply OpenDX(formerly IBM Data Explorer)NCL(NCAR Command Language)ncBrowse(netCDF File Browser)HDF(Hierarchical Data Format)interface GrADS(Grid Analysis and Display System)FERRET NetCDF Toolbox for MATLAB-5 IDV(Integrated Data Viewer)Originsurfer,常用绘图软件，及GrADS优缺点,matlab画图和grads画图的比较 总体上来说matlab画图比grads好 matlab的优点：（1），matlab是用对象属性来控制图形中的每一个部分，因此在对图片进行操作是相当容易。（2），它画二维曲线图时，画的图非常光滑，视觉效果很好。（3），matlab画二维等值线图时，投影方式相当丰富，基本覆盖各行业所需。它的底图（地图）比较漂亮，而且还可以用海拔高度函数是不同海拔用不同的颜色表示出来。（4），它可以画三维图，图形非常漂亮。grads在这点目前是办不到的.Matlab的缺点：绘图速度较慢，当数据文件较大时，读入很容易内存溢出远程处理比较麻烦。grads的优点;（1），最多能读入20个数据文件，数据读取迅速，出图速度快，但是要写描述文件和gs文件比较麻烦。（2），在画二维等值线图时，图形非常平滑，这点是我到现在还没大规模用matlab的主要原因。（3），它包含了很多气象上所需要的函数，在gs文件里可以直接调用，因此在用它时也省了不少事。,Grads在win/Linux下的下载、安装、运行,程序：Windows：Linux：环境变量：Windows：新建变量GASDIR c:PcGrADSdatGASCRP c:PcGrADSlib修改PATH变量,增加路径c:PcGrADSwin32;c:PcGrADSlib;c:PcGrADSdat;Linux：在.cshrc 文件中添加:#执行文件路径：setenv PATH/usr/local/grads/bin:$PATH#字体和底图文件路径：setenv GADDIR/usr/local/grads/dat#脚本库的路径：setenv GASCRP/usr/local/grads/lib,Grads网络资源,GrADS主页和中文讨论区 或(软件下载/联机文档/email讨论区)(中文手册下载/中文讨论区/脚本和其它资源)GrADS实用文档 LASG编GrADS实用手册中文,较全面,可用来自学 GrADS reference card英文,绘图命令和函数速查 GrADS scripting language reference card英文,脚本语言语法速查 GrADS Documentation Index联机英文索引,最全面,运行GrADS,打开数据：ga-open*.ctlga-sdfopen*.ncga-xdfopen test.ddf查询信息：queryga-query ctlinfoga-q dimsga-q file关闭数据文件ga-close 3ga-close 2Ga-close 1设置维数环境ga-set x 1 73 ga-set lat-20 60 ga-set lev 850 ga-set time 00z01Apr2007 26LASG动力论坛版三.运行GrADS绘图软件包 设置维数环境ga-set x 1 73 ga-set lat-20 60 ga-set lev 850 ga-set time 00z01Apr2007,显示图形ga-display u ga-d u;v.2 ga-d u;v.2;mag(u,v.2)保存图形ga-enable print a.gmf ga-print ga-disable print 清屏ga-clear ga-c 重置缺省环境,不关闭数据文件ga-reset 重新初始化,关闭reinit,命令顺序grads b|l|p|copen/runsetdisplaydrawprintquit/reinit,显示类型,bar 直方图line 折线图errbar 标误差线linefill 两线之间填色scatter 散点图grfill 按网格填色(马赛克图)fgrid 按set fgvals指定值填网格颜色grid 网格填值findstn 寻找最近的站点model 站点填图stnmark 标记站点位置value 标站点值contour 等值线图shaded 填色等值线图vector 风矢量图barb 风羽stream 流线图stat 输出统计数据print 十进制屏幕输出(set prnopts)fwrite 输出二进制文件(set fwrite)tserbarb 单站的风羽时间序列tserwx 单站天气现象时间序列wxsym 在站点上标记天气符号(set wxopt),坐标设置,title/xlab/ylab 和坐标轴/刻度/标值-画图形区的外框-X和Y轴互换-Z轴取对数坐标(气压不等距)-(1.9)纬度1-D变化时纬度取cosine坐标-省略时间刻度的年份/年月-绘图区内部网格线属性-控制set gxout grid的内部网格线属性-刻度间隔-刻度及标值的起/止位置及间隔-Y轴变量取值范围(或X轴)-坐标轴相对于绘图区的位置-坐标刻度和标值的属性-等值线标值的属性-列举坐标刻度-列举坐标标值-X轴反向(或Y轴)-坐标刻度的标值方式-写坐标标题-图的标题-,set annot#color(1)thick(6)set frame on|off|circleset xyrev onset zlog on|off|swap|undefine set coslat on|off set tlsupp year|month set grid on|off|horizontal|vertical linestyle#colorset gridln auto|off|#color set xlint intset xaxis starrt end intset vrange vlo vhiset xlpos offset(0)b|t set xlopts#color(1)thick(4)size(0.12)set clopts#col(-1)thick(-1)size(0.09)set xlevs x1 x2 x3set xlabs b1|b2|b3|set xflip on|off set xlab on|off|autdraw title title_string,set ylintset yaxisset vrange2set ylposl|rset yloptsset ylevsset ylabsset ylabset yflipdraw ylab,自由图形、符号,设置线属性-画线-画符号-画空心|实心方框-画多边形-画多边形-画天气符号-设置字体-写标题-设置字符串属性-设置字符串大小-写字符串-,set line#color#style#thickdraw line x1 y1 x2 y2draw mark marktype x y sizedraw rec|recf xlo ylo xhi yhidraw ployf x1 y1 x2 y2 xn yndraw mappoly lon1 lat1 lon2 lat2 lon1 lat1draw wxsym#symbol x y size#color#thickset font 050draw title title_stringset string#color tl|tc|tr|l|c|r|bl|bc|br set strsiz width heightdraw string x y string,天气符号示例,draw wxsym#symbol x y size#color#thick,q w2xy 135 20 ga-X=4.25 Y=4.91667draw wxsym 40 5 4 0.5 2 2,地图投影、地图数据集,地图投映方式set mproj latlon|scaled|nps|sps|robinson|orthogr|mollweide|lambert|off 地图数据名set mpdset lowres|mres|hires|地图属性设置set map auto|#color#style#thickness 地图绘制开关set mpdraw on|off 行政区界开关set poli on|off 投映参考经纬度set mpvals off|lonmin lonmax latmin latmax,地图投影方式：set mproj latlon|scaled|nps|sps|robinson|orthogr|mollweide|lambert|off,-nps-sps,-orthogr,mollweide-,robinson-,地图数据集,图形、数据输出,图形输出,1,可用于批处理文件中enable print filename.gmfprintdisable printprintim png|gif filename xNNN yNNN white|black 2，不能用于批处理文件中wi filename.fmt.fmt 可为任意图形格式后缀注意图形窗口不能被遮挡！outxwd filename.xwd,数据输出,1，输出成二进制文件set gxout fwriteset fwrite filenamed vardisable fwrite2，输出为ASCII文件file=output.txtset gxout printd temrc=write(file,result)rc=close(file),数据分析与诊断函数的使用,数学函数绝对值abs(expr)余弦/反余弦cos(expr)acos(expr)正弦/反正弦sin(expr)asin(expr)正切/反正切tan(expr)atan2(expr1,expr2)expr1/expr2指数exp(expr)对数log(expr)log10(expr)幂函数pow(expr1,expr2)expr1expr2全风速mag(u_expr,v_expr)平方根sqrt(expr)统计函数：网格处理函数 滤波平滑函数 差分函数 物理量计算函数 站点数据专用函数,统计函数：,加权/不加权平均-加权/不加权区域平均-剔除值后时间求平均-加权/不加权区域和-加权/不加权区域和-质量加权垂直积分-普通积分-最大/小值-最大/小值所在格点-空间相关/回归-时间相关/回归-,ave(expr,dim_expr1,dim_expr2,tincr,|-b)mean()aave(expr,xdim1,xdim2,ydim1,ydim2|global|g)tmave(undefine_or_weight_expr,expr,texpr1,texpr2)sum(expr,dim_expr1,dim_expr2,tincr,|-b)asum(expr,xdim1,xdim2,ydim1,ydim2|global|g)vint(ps_expr,expr,top)gint(expr)Max(expr,dim1,dim2,tincr)min()Maxloc(expr,dim1,dim2,tincr)minloc()scorr(expr1,expr2,xdim1,xdim2,ydim1,ydim2|global|g)tcorr(expr1,expr2,tdim1,tdim2),网格处理函数,常数 const(expr,const,|-u|-a)剔除 maskout(expr,mask_expr)跳点skip(expr,skip_x,skip_y)加权9点平滑 smth9(expr)1-D 3点平滑中央差分 cdiff(expr,dim),物理量计算函数,虚温,相对湿度求比湿-虚温,相对湿度求温度-垂直涡度-中央差求水平散度-全风速-,tvrh2q(tv_expr,rh_expr)gtvrh2t(tv_expr,rh_expr)hcurl(u_expr,v_expr)hdivg(expr1,expr2)mag(u_expr,v_expr),命令批处理文件：Grads Script:*.gs,注释语句:*赋值:交互命令：say/prompt/pull条件控制：if/else/endif循环语句：while endwhile/break/continue,*this is a sample script open my_sst_dataset.ctl set lat-30 30 set lon 90 300 display sstsay She said it is.Prompt Enter min and max longitudes:pull minlon maxlon if expressionscript record.elsescript record.Endif while expressionscript record.endwhile,GrADS自带脚本库,加注色标/图例cbar.gs 条状色标cbarn.gs 两头尖的条状色标cbarc.gs 扇形色标cbar_l.gs 折线图图例cbar_line.gs 折线图图例地图底图basemap.gs 海洋/陆地填色makebg.gs 生成具有地形纹理的底图map.gs 预定义多个可选的地图范围颜色定义cmap.gsdefine_colors.gs,插值isen.gs 由等压面插值到特定等熵面pinterp.gs 由eta/sigma/等压面等插值到指定等压面zinterp.gs 由eta/sigma/等压面等插值到指定等高面物理量计算等mconv.gs 计算水汽辐合plotskew.gs 绘制T-logP图sweat_index.gs 用相对湿度,温度和风计算凝结指数traj.gs 在水平面上绘制当前时刻向前/向后的轨迹线格式转换lats4d.gs 输出NetCDF,HDF-SDS 或GRIB 格式数据动画/缩放xanim.gs 动画控制zoom.gs 图形放大,GrADS数据数据类型与格式,支持的数据格式：netCDF数据格式HDF 数据格式Binary数据格式站点数据格式GRIB数据格式,netCDF数据格式,运行gradsncga-sdfopen path/filename.nc,无格式的Binary数据格式,自西向东 从南到北 从低层到高层-变量1到变量n-时间从小到大这是GrADS所用的数据排列方式,do it=1,mtwrite(10,rec=it)(var(I,j,k,it),i=1,mx),j=1,my),k=1,ml)enddo,parameter(mx=144,my=73,ml=7,mt=12)dimension slp(mx,my,mt),Dimension uwnd(mx,my,ml,mt),vwnd(mx,my,ml,mt)open(10,file=sample.grd,form=unformatted,access=direct,recl=mx*my)irec=1do it=1,mtwrite(10,rec=irec)(slp(i,j,it),i=1,mx),j=1,my)irec=irec+1do ilev=1,mlwrite(10,rec=irec)(uwnd(I,j,ilev,it),i=1,mx),j=1,my)irec=irec+1enddodo ilev=1,mlwrite(10,rec=irec)(vwnd(I,j,ilev,it),i=1,mx),j=1,my)irec=irec+1enddoenddoclose(10)end,DSET sample.datTITLE Sample Model DataDTYPE station UNDEF 0.10000E+16XDEF 144 linear 0 2.5YDEF 73 linear-90 2.0ZDEF 7 levels 1000 900 800 700 500 300 100TDEF 12 linear 00z01apr1985 6hrVARS 3slp 7 99 sea level pressure u 7 99 U component v 7 99 V componentENDVARS,模式数据转化为nc数据,MM5toGrADS.tar.gz:WRF2GrADS:,站点数据的准备以及站点图形的绘制,set gxout 图形类型findstn 寻找最近的站点model 站点填图stnmark 标记站点位置value 标站点值tserbarb 单站的风羽时间序列tserwx 单站天气现象时间序列wxsym 在站点上标记天气符号(set wxopt),台站型离散资料画图的一般原则,1、把台站资料r.dat写成二进制文件r.grd2、编写r.grd数据的数据说明文件r.ctl3、用r.ctl生成r.map文件4、生成插值用的背景格点文件5、编写grid.grd的.ctl文件：grid.ctl6、把r.grd插值到grid.grd的网格点上，并画图,第一、二步：把台站资料*.dat写成二进制文件*.grd，并编写相应ctl文件,Fortran 程序片段,对应的数据描述文件,dset sample.grdtitle Sample Model Datadtype station stnmap sample.mapundef 0.10000E+16xdef 144 linear 0 2.5ydef 73 linear-90 2.0zdef 7 levels 1000 900 800 700 500 300 100tdef 12 linear 00z01apr1985 6hrvars 3slp 7 99 sea level pressure u 7 99 U component v 7 99 V componentendvars,parameter(mx=144,my=73,ml=7,mt=12)dimension slp(mx,my,mt),Dimension uwnd(mx,my,ml,mt),vwnd(mx,my,ml,mt)open(10,file=sample.grd,form=unformatted,access=direct,recl=mx*my)irec=1do it=1,mtwrite(10,rec=irec)(slp(i,j,it),i=1,mx),j=1,my)irec=irec+1do ilev=1,mlwrite(10,rec=irec)(uwnd(I,j,ilev,it),i=1,mx),j=1,my)irec=irec+1enddodo ilev=1,mlwrite(10,rec=irec)(vwnd(I,j,ilev,it),i=1,mx),j=1,my)irec=irec+1enddoenddoclose(10)end,第三步：用stnmap 外部命令生成.map文件,方法1、ga-!stnmap enter stn ctl filename:e:/data/sample.ctl方法2、!stnmap i e:/data/sample.ctl 特别注意：GrADS中的路径表达方法:1,引用外部命令时，用unix表达法反斜杠/，例如用!stnmap i e:/data/sample.ctl，否则要用!stnmap i e:datasample.ctl2,在grads的命令提示符下输入的内部命令，如open等，和ctl描述文件中的路径则不区分/,和,例如：ga-open e:datasample.ctl 等同 ga-open e:/data/sample.ctldset e:datasample.grd 等同 dset e:/data/sample.grd,第四步、生成插值用的背景格点文件,因为GrADS只识别格点数据，所以需要把站点上的资料通过插值函数插值到某个格点文件上才能被识别。根据r.dat提供的信息来写格点文件grid.grd（二进制）要求：（1）grid.grd文件的精度要高于或等于r.dat的精度（2）grid.grd文件的范围要大于或等于r.dat的范围（3）grid.grd文件的每个点上均赋值1（范围与r.dat相同）（4）当grid.grd文件的范围大于r.dat的范围时，除了相同范围的每个点赋值为1外，其余外围的点也都赋值为0，（也可以赋值为1，不过这样做最后在站点数据插值到网格点时，会插值到所期望的网格区域外）（5）插值网格文件的ctl文件的时间一定要与数据资料的时间一致！否则会出现Warning from OACRES:Less than two stations的错误信息提示,dset e:datachina.grdtitle gridundef 999.9xdef 620 linear 73.0 0.1ydef 359 linear 18.0 0.1zdef 1 levels 1000 tdef 1 linear jan1951 1yrvars 1g 0 99 grid data1endvars,相应的ctl描述文件：,第五步、将站点资料插值到预先定义的网格上，并绘图,open e:datachina.ctlopen e:datasample.ctlset lon 70 140set lat 15 55define aa=oacres(g.1,t1.2,50,40,30,20,10,5,2)define aaa=maskout(aa,g.1-0.5)set gxout shadedd aaaset gxout contourd aaaprintim e:datasample.png white,