《GrADS处理台站数据精讲》由会员分享，可在线阅读，更多相关《GrADS处理台站数据精讲（6页珍藏版）》请在金锄头文库上搜索。
1、Nanjing university of information science& technology 阿木 废话不说，开始讲解第一步：处理站台id文件： 文件处理成三列数据，第一列为站台号，第二列为纬度，第三列为经度 注意，如：站号 纬度 经度50136 52.97 122.5250246 52.35 124.7250349 51.70 124.3350353 51.72 126.65上述数据中50136连同后面的空格要保证一共有8位（这点不容质疑）第二步：处理数据 一般要将所要的要素，如日照时间，单个要素放到一个txt文件中，数据的顺序是：同一个时间点371个站数据，然后放下一个时间点
2、的371个要素值，如此直到最后时刻。第三步：处理成对应格式数据grd文件案例程序（371个站18263天的数据）：说明：这里的rct_50years_371_5vars50year_371id.txt 为第一步生成的台站id文件 rct_50years_371_5vars371_sd_1956_2005.txt 为第二步生成的数据文件注意：1、数据必须设成实型的，站号必须为8位字符型（与第一步照应）； 2、数据存放依然是先站，后时次，即最外层为时间循环 3、每一行数据：站号纬度经度相对时次记录标志nflag要素数据 4、程序中tim为相对时次，指的是相对于一天的世界时，如果数据本身是日数据，那
3、么time就取0.0，一般都是这样，除非是有08时、20时等多个时刻，才考虑time不为0.0，此处设为0.0恒定，元数据为日数据 5、程序中nlev标志一个时刻记录开始，结束时赋值为0，对应程序； 6、nflag为有没有地面报，有为1，否则为0，此处为地面观测相对湿度，为1 7、循环中每个时次的371个站数据输完后，要加一行标志记录（红色） 8、生成的文件为二进制grd文件，所以为write（30），没有格式说明；直接加文件号 9、第一二三步数据的顺序一定要与站号顺序对应。program ex00implicit noneinteger nlev,nflag,i,jreal lat(371)
4、,lon(371),tim,sd(371)character*8 stid(371)tim=0.0nflag=1open(10,file=D:Program FilesFORTRAN6.0txt_gridrct_50years_371_5vars50year_371id.txt)open(20,file=D:Program FilesFORTRAN6.0txt_gridrct_50years_371_5vars371_sd_1956_2005.txt)open(30,file=D:Program FilesFORTRAN6.0txt_gridrct_50years_371_5vars371_s
5、d_1956_2005.grd,form=binary)do i=1,371 read(10,*) stid(i),lat(i),lon(i) enddoclose(10)do i=1,18263 nlev=1 do j=1,371 read(20,*) sd(j) write(30) stid(j),lat(j),lon(j),tim,nlev,nflag,sd(j) enddo nlev=0 write(30) stid(371),lat(371),lon(371),tim,nlev,nflag write(*,*) Finish writting:,ienddoclose(20)clos
6、e(30)end第四步：写相应数据grd的ctl，如下dset c:371yue371_sd_1956_2005.grddtype stationstnmap c:371yuesd50.mapundef 32766title sun hourtdef 18263 linear 01jan1956 1dyvars 1sd 0 99 sunhour(0.1h)endvars第五步：写标准的格点文件用于后期插值这里经度为1*1，70E-140E,15N-55-N注意：grd的经纬度范围要大于格点的范围，时间可对应。parameter(nx=71,ny=41,nt=600)real lat(ny),l
7、on(nx)real s(nx,ny)open(1,file=grid.grd,form=binary)lat(1)=15.0lon(1)=70.0do j=1,ny-1 lat(j+1)=lat(j)+1.0enddodo i=1,nx-1 lon(i+1)=lon(i)+1.0enddodo k=1,18263 do i=1,nx do j=1,ny s(i,j)=1enddo enddoenddowrite(1) send第六步：写出对应插值grd的ctldset c:371yue371_sd_1956_2005.grddtype stationstnmap c:371yuesd50.m
8、apundef 32766title sun hourtdef 18263 linear 01jan1956 1dyvars 1sd 0 99 sunhour(0.1h)endvars第七步：生成数据地图文件台站的数据文件，需要地图文件才能读取方法： 在grads启动显示图形窗口后，输入 ! stnmap 回车 输入：c:371yue371_sd_1956_2005.ctl 回车 然后就在对应路径下生成相应的地图文件了，接下来就可以用grads画台站图了。如下，台站旁有数据（不是0）才说明你前面的处理对了由于grads能画的台站图比较局限，有下面的差值成格点数据第八步：插值，画图案例：gs文件
9、：reinitopen c:371yuegrid.ctlopen c:371yueSD50.ctlset annot 5 8set xlopts 3 5 0.16set ylopts 3 5 0.16set parea 0.6 10.5 0.9 8.0set grads offset grid offset lat 15 55set lon 70 140set t 1 18263define a=oacres(g(t=1),SD.2)define aa=smth9(a)set t 1save=ave(aa,t+0,t=18263)set t 1set gxout shadedset grads
10、 offcnbasemap savecbarndraw title 50nianave_SD.pngprintim c:50nian_SD.png white说明：1、这里主要插值在红色处，g为第五步生成的grd文件的变量，SD为第三步生成grd文件变量，oacres（g（t=1），SD.2）为将2号文件的SD量插值到1号文件的g上，插值半径没写，取为默认，笔者认为插值半径默认就好，如果自己随便调的话，容易出其他问题； 2、这里画图命令是用cnbasemap save，如果用原来的 d save，画出的图会超出中国边界，这个命令可以自动将边界外线抹除，图很漂亮，但是，要用这个命令，需要装1.9版本的grads，还要装补丁。 3、有需要补丁的朋友或者对文章有什么疑问的发邮件给我,邮箱：附图二张： 1
