grads处理多个ctl文件及nc文件.docx

grads处理多个ctl文件和nc文件  
2011-10-10 21:03:59|  分类： grads学习 |  标签： |举报 |字号大中小 订阅

  下载LOFTER 我的照片书  |
用grads处理多个相同格式的数据时若单个单个处理非常麻烦，当文件非常多的时候是单个处理是不实际的。下面介绍一种方法;
   第一步，在这种情况下可以重新写一个ctl描述文件，其文件变量都和已知的ctl相同，若原来的n文件只是时间不同，那么新描述文件的时间维数是所有原文件的时间的和。同样，若其他维数不同时也用同样的方法处理。 
   第二步，在第一行之后添加一行 options template 表示多个时间序列原始数据文件想用一个描述文件统一地描述。这些原数据的原文件名由dset定义的形势命名文件名。
   第三步，修改dset 的文件名。原路径不变，把文件名用％表示。其中：
   ％y2  代表两位数年
   ％y4  代表四位数年
   ％m1  代表一位或者两位数的月
   ％m2  代表两位数月（用0补齐1位数）
   ％mc  3个字符月份的缩写
   ％d1  1或2位天
   ％d2  两位天
   ％h1  1或者2位时
   ％h2  2位时
例如：
   原文件其中之一的文件名为gdas2006050812f00，且所有文件只有天和时的变化
   那么新描述文件的文件名为：gdas200605%d2%h2f00

   另外如果源文件里有index项的话，需要修改其idx的文件名，假设改成fnl.idx。并用在dos下用gribmap函数生成一个新的idx文件。gribmap -e -i fnl.ctl（加绝对路径）
open fnl.ctl就可以打开所有文件。

**********************************************************************************************************
若想要提取从1951-2006年56年nc文件中的某些数据，一个一个处理非常麻烦，这里介绍种较为简易的方法。例如想提取6-8月的位势高度资料。
'reinit'
t5=1951
*作文件名循环
while(t5<=2006)
 'set gxout fwrite'
 'set fwrite D:\sichuan\hgt1\'%t5%'.dat'
 'sdfopen e:\ncep1\hgt\hgt.'%t5%'.nc'
 t3=t5-1950
*判断是否为闰年
 if（t3=2|t3=6|t3=10|t3=14|t3=18t3=22|t3=26|t3=30|t3=34|t3=38|t3=42|t3=46|t3=50|t3=54)
  to=153
 else
  to=152
 endif

   t4=to+91
   while(to<=t4)
     'set t 'to
      t1=1
      while(t1<=12)
        'set z 't1
        'set lon 80 140'
        'set lat 15 55'
        'd hgt'
         t1=t1+1   
      endwhile

      to=to+1
   endwhile
*这里必须先观点上述运行的文件，grads最多同时可以打开20个文件左右。
   'reinit'
 t5=t5+1
endwhile
'reinit'
这样可以提取你想要的年数据，然后你大可运用fortran对数据进行随心所欲的处理。

能否直接生成一个文件还正在探索中。

批量读取nc数据，用你的方法成功了，谢谢！！！直接配个批量描述的ctl就可以了
有一批nc数据，一个月一个文件，现将文件名改为：197901.nc，197902.nc，依次类推，对二进制的数据知道写ctl文件来进行批处理运算，那么nc数据应该怎么做呢？试过了写ctl文件，sdfopen ***\%y4%m2.nc，
year=1978
while(year<=2011)
month=01
  while(month<=12)
'sdfopen ***\'year''month'.nc'
...
    month=month+1
  endwhile
year=year+1
endwhile

实我也是糊里糊涂的解决了。。。ctl文件如下：
dset  ^%y4%m2.nc
undef 1e+15
options template
title MERRA data
dtype NetCDF
ydef  144 linear -90 1.25
xdef  288 linear -180 1.25
zdef   21 levels 1000 975 950 925 900 875 850 825 800 775 750 725 700 650 600 550 500 450 400 350 300
tdef  396 linear 00Z01JAN1979 1mo
vars 3
qv 21 t,z,y,x Specific humidity
u 21 t,z,y,x Eastward wind component
v 21 t,z,y,x Northward wind component
endvars
然后open ***.ctl就行了，之前的问题是打不开ctl文件，怎么改也不行，后来换了台机子就好了。。。所以我说我也不知道怎么回事，希望对你有帮助。

我之前就这样做的，能打开ctl文件，但是d之后，都显示all undefined values，我的ctl如下，麻烦帮我看看哪里错了？
dset  ^%y4%m2.nc
undef 1e+15
options template
title MERRA data
dtype netcdf
xdef  288 linear -179.375 1.25
ydef  144 linear -89.375 1.25
zdef   21 levels 1000 975 950 925 900 875 850 825 800 775 750 725 700 650 600 550 500 450 400 350 300
tdef  396 linear 00Z01JAN1979 1mn
vars 3
qv 21 t,z,y,x Specific humidity
u 21 t,z,y,x Eastward wind component
v 21 t,z,y,x Northward wind component
endvars
可以的，但文件名一定要连续
这个时间长度一定要和你的文件对应好

最近发了一个利用grads批量合并nc文件的帖子  很多人也都问怎么批量描述，，很多人也都去自己尝试，尝试过程中出了各种问题。推荐兰溪给nc文件写ctl的帖子
“SDF file has no discernable X coordinate”的同学们，非常推荐这个帖子。肯定很多人已经看过这个帖子，也照着做过，还是出了问题。

其实论坛还有很多其他有关ctl批量描述nc的帖子，大家参考一下，对照自己出现的问题，应该大部分情况下都能解决。但是有一些人很着急，连一个nc文件的描述文件都写不对，就直接批量的，那肯定只会出更多的错误，还一时半会儿改不对。只会造成时间的浪费。
  我不是来说教的，只是我觉得什么事儿都是循序渐进的，不要那么浮躁。俗话说磨刀不误砍柴工，其实真的是这么回事。就说批量描述nc文件，它也是在正确描述一个nc文件的基础上来修改的。只要描述一个的出来了，往上加一两个语句，用正确的替代格式替代了文件名，再改一下时间长度就出来了。那能正确描述一个nc文件的ctl就是你磨得很锋利的砍柴刀，有了它后面的柴就相当好砍了。
  下面我就以ncep逐日再分析资料为例说一下nc文件的描述和批量描述。文件名连续，hgt.1948.nc
hgt.1949.nc  hgt.1950.nc`````````
首先看一下它自带的ctl，如下图红色圈起来的部分，就不详细解释了，对grads有一定了解的人，看一下就知道什么
自带ctl.PNG (29.64 KB, 下载次数: 9)
下载附件  保存到相册
2013-5-22 22:42 上传

下面就可以根据这个自己编写一个ctl文件来描述这个nc文件了。
  照着自带的写下来先，存成1948.ctl然后使用open命令打开，画图，发现出来的图和缺测值设置错误时差不多。那就想方设法的改缺测值，比如用在grads中使用q attr查看有一句missing_value 32766，修改缺测值，再画图，还是那样；q undef出来的是-999000000.000000，再改，再画还是那样。无论怎么改，都还是那样。
  如果你这么反复折腾，最后还没发现问题，那就说明你没有好好利用论坛的搜索功能，也没有看到兰溪的帖子，也没看到黎大叶子的用Fortran批量为nc写ctl的帖子

  但是不用着急，我看见了，其中最重要的是打开netcdf格式数据的描述文件是需要用xdfopen命令的，那就先要去看看xdfopen能打开的ctl需要怎么写http://grads.iges.org/grads/gadoc/gadocindex.html。
  看完了，差不多就能明白了，有这么多前人的成果那就照着修改呗，最终我修改出来了，图画的也很正常了。写出来的如下
1. dset F:/ncep/daily/hgt.1948.nc
2. title mean daily NMC Reanalysis （1948）
3. undef -999
4. xdef lon 144 linear 0 2.5
5. ydef lat 73 linear -90 2.5
6. zdef level 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
7. tdef time 366 linear 00Z01Jan1948 1440mn
8. vars 1
9. hgt=>hgt 17  t,z,y,x  mean Daily Geopotential height
10. endvars
复制代码
  光看着正常不行啊，需要和原始的图对比验证了才能确定是对的吧。所以在grads里面用sdfopen命令打开hgt.1948.nc画第一层第一个时次的图，再用xdfopen打开你编写的ctl，也画第一层第一个时次的图看看。我擦！！！！对不上啊，看来还是有问题啊。说明上面的ctl是有问题的，还得改进才行。  
  其实到这步已经基本成功了，仔细看看叠加起来的那两张图，几乎是一样的，只是南北的方向是反的。那就好办了，在ctl里加上options yrev，告诉grads南北要反向不就行了么。于是最终的ctl出来了，如下
1. dset F:/ncep/daily/hgt.1948.nc
2. title mean daily NMC Reanalysis （1948）
3. options yrev
4. * yrev表示y轴反向
5. undef -999
6. xdef lon 144 linear 0 2.5
7. ydef lat 73 linear -90 2.5
8. zdef level 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
9. tdef time 366 linear 00Z01Jan1948 1440mn
10. vars 1
11. hgt=>hgt 17  t,z,y,x  mean Daily Geopotential height
12. endvars
复制代码

再用xdfopen打开画图，这回就一模一样了啊。 说明成功了！

那下面的批量描述就太简单了，比如我要批量描述1948和1949两年的，算一下一个闰年一个平年，一共有时次366+365=731，那么就修改ctl吧
dset F:/ncep/daily/hgt.%y4.nc
title mean daily NMC Reanalysis
options template
options yrev
* yrev表示y轴反向
*undef -999
xdef lon 144 linear 0 2.5
ydef lat 73 linear -90 2.5
zdef level 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
tdef time 731 linear 00Z01Jan1948 1440mn
vars 1
hgt=>hgt 17  t,z,y,x  mean Daily Geopotential height
endvars
就这么简单，只要用%y4代表四位年，把总的时次改成731，再增加一句options template就可以了。然后就可以利用xdfopen命令打开画图和原来的对比一下，一样一样的吧，可以批量描述了吧！
注意：有人可能注意到我把缺测那一项给注释掉了，其实这个是完全不需要的，去掉或者改成任何值，都不影响。这是问什么？不细说了，因为我没仔细研究（比较晚了，我要去睡觉了，以后再说，或者有人知道可以告诉我）

其他的资料都同理了，思路都差不多，变通一下就可以了。
1. 首先编写可以正确描述一个资料的ctl，并能正确出图。
2. 修改该ctl里dest后的文件名，使用合适的替代格式替代（具体的格式看实用手册）
3. 增加options template，这个是批量描述必须加上的一句
4. 计算出所要批量描述的文件的总时次，修改ctl里的总时次
5. 画图进行验证，如果不对，再根据具体情况做出相应的修改
6. 大功告成，可以用了
7. 需注意的一点就是options yrev，这是可选项，要根据资料实际情况来使用。

今天这个帖子有点儿罗嗦了，大家捡有用的看就行。其实说那么多，无非是想和一些人说我们这个年代的人已经是站在前人的肩膀上了，很多东西别人都已经有了经验，分享给你你就要好好利用，不要把时间都浪费在发帖和等回复上面。很多东西，只要自己有一定的基础，在一定合理的范围内就会找到那个答案的·····
[GrADS] grads批量处理nc文件求区域平均
[复制链接] |关注本帖
· 取消最新回复
· 取消置顶回复
· 取消最新编辑

stefan_scofield
stefan_scofield 当前离线
积分
4181
贡献
精华
在线时间
小时
注册时间
2013-7-18
最后登录
1970-1-1
窥视卡 雷达卡

电梯直达
楼主
海温资料是1965-2013年的逐月资料（每年的4,5,6月），我需要求每年4,5,6月nino3.4区域的平均海温，然后三个月再平均，就是每年得一个值。

win7

批量处理nc文件，合成一个dat文件时遇到问题，请大家帮指点一下！
'reinit'
'set gxout fwrite'
'set fwrite g:\mls\1.20.32.dat'
t5=1
*作文件名循环
while(t5<=20)
'sdfopen g:\mls\h2o.'t5'.nc'
         'set x 1 144'
         'set y 1 91'
         'set t 1'
         z=1
         while(z<=25)
         'set z 'z
          'd h2o.'t5
          z=z+1
       endwhile
t5=t5+1
endwhile
'disable fwrite'

grads提取某个范围内的资料时注意的一个问题
热度 4已有 524 次阅读2011-11-19 16:10 | 资料
这两天被一个资料的提取程序弄的晕乎乎的，今天终于发现问题所在了。
提取资料的时候，经度的范围即使原来就是全球的，也不能写成 set lon 0 360，得写成set lon 0 357.5，就这个2.5（格点间距）把我折腾了这么久，下次一定要注意了！

海温资料是1965-2013年的逐月资料（每年的4,5,6月），我需要求每年4,5,6月nino3.4区域的平均海温，然后三个月再平均，就是每年得一个值。报错：Data request warning: request is completely outside file limits

gs如下：
'reinit'
'set gxout fwrite'
'set fwrite d:/sst/4-6nino.grd'
yy=1965
while(yy<=2013)
'sdfopen d:/sst/ersst.'yy'04.nc'
'sdfopen d:/sst/ersst.'yy'05.nc'
'sdfopen d:/sst/ersst.'yy'06.nc'
;'set z 1'
;'set t 1'
;'set lat -5 5'
;'set lon 170 240'
'a=aave(sst.1,lon=170,lon=240,lat=-5,lat=5)'
'b=aave(sst.2,lon=170,lon=240,lat=-5,lat=5)'
'c=aave(sst.3,lon=170,lon=240,lat=-5,lat=5)'
's=(a+b+c)/3.0'
'd s'
'close 3'
'close 2'
'close 1'
yy=yy+1
endwhile
'disable fwrite'
;

我编写的批量描述文件如下：
DSET d:/TRMMdata/3b42/9801/3b42.2000%m201.00.7a.nc
TITLE TRMM data of year 2000
DTYPE netcdf
OPTIONS template
UNDEF -9.99E33
XDEF 144 LINEAR 0 2.5
YDEF 20 LINEAR 0 2.5
ZDEF 1 LEVELS 1000
TDEF 12 LINEAR 01JAN2000 1mo
VARS 1
pcp=>pcp 0 t,z,y,x precipitation
ENDVARS

gs文件如下：
'reinit'
'open d:/trmmdata/work/2000.ctl'
'set lat 0 60'
'set lon 0 30'
'set lev 1000'
'd ave(pcp,t=1,t=12)'
'reinit'

 用GrADS转换nc数据
                                                      ——by Artmunich from
    kittyhare 在一篇帖子里向我们初步介绍了使用GrADS转化nc为grd格式（摸我去看看~），我想能不能再详细点，于是借鉴以前自己见到的一个教程，更详细的介绍一下这方面的知识，希望大家指正。
   先给出一个单层的二进制文件转化的gs

'reinit'
'sdfopen f:/data/nc/1980/air.1980.nc'
'set gxout fwrite'
'set fwrite f:/data//air.1980.bin'
'set lon 0 357.5'
'set lat -90 90'
'set lev 1000'
'set t 1 640'
'd air'
'reinit'

    但我们知道nc文件是按照经度、纬度、高度、变量、时次顺序排列，要转nc文件，需要靠循环，由于高度的不连续性，我们可以在时间循环里面把每层的高度写出来。
'reinit'
'sdfopen f:/data/nc/1980/air.1980.nc'
'set gxout fwrite'
'set fwrite f:/data//air.1980.bin'
'set lon 0 357.5'
'set lat -90 90'
t=1
while(t<=365)
'set t' t
'set lev 1000'
'd air'
'set lev 925'
'd air'
*You can continue to write hgt
t=t+1
endwhile
'reinit'
   当然，如果用z坐标系，这样z就是从1到17的，在时间的循环里嵌套z=1;while(z<=17)我认为也是可行的，不过自己没试过，有用过的朋友可以告诉一下结果。
高度循环的：
while(t<=8)
'set t 't''

z=1
while(z<=21)
'set z ' z''
'set lon 90 125'
'set lat 10 35'
'define t=tmpprs'
'define rh=rhprs'
'define es=(6.112*exp(17.67*(t-273.15)/(t-29.65)))'
'define q=rh*(0.62197*es/(lev-es))/100.'
'define e=lev*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(t)-log(e)-4.805)'
'define theta=t*pow((1000./lev),(0.2854*(1.0-0.28*q)))'
'define thetse=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
'd thetse'
z=z+1
endwhile
t=t+1
endwhile