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

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

《grads处理多个ctl文件和nc文件.docx》由会员分享，可在线阅读，更多相关《grads处理多个ctl文件和nc文件.docx（16页珍藏版）》请在冰豆网上搜索。

grads处理多个ctl文件和nc文件

grads处理多个ctl文件和nc文件  

2011-10-1021:

03:

59|  分类：

grads学习|  标签：

|举报|字号大中小 订阅

  下载LOFTER我的照片书  |

用grads处理多个相同格式的数据时若单个单个处理非常麻烦，当文件非常多的时候是单个处理是不实际的。

下面介绍一种方法;

  第一步，在这种情况下可以重新写一个ctl描述文件，其文件变量都和已知的ctl相同，若原来的n文件只是时间不同，那么新描述文件的时间维数是所有原文件的时间的和。

同样，若其他维数不同时也用同样的方法处理。

  第二步，在第一行之后添加一行：

optionstemplate表示多个时间序列原始数据文件想用一个描述文件统一地描述。

这些原数据的原文件名由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-ifnl.ctl（加绝对路径）

openfnl.ctl就可以打开所有文件。

**********************************************************************************************************************************

若想要提取从1951-2006年56年nc文件中的某些数据，一个一个处理非常麻烦，这里介绍种较为简易的方法。

例如想提取6-8月的位势高度资料。

'reinit'

t5=1951

*作文件名循环

while（t5<=2006）

 'setgxoutfwrite'

 'setfwriteD:

\sichuan\hgt1\'%t5%'.dat'

 'sdfopene:

\ncep1\hgt\hgt.'%t5%'.nc'

 t3=t5-1950

*判断是否为闰年

 if（t3=2|t3=6|t3=10|t3=14|t3=18|t3=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）

    'sett'to

     t1=1

     while（t1<=12）

       'setz't1

       'setlon80140'

       'setlat1555'

       'dhgt'

        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

undef1e+15

optionstemplate

titleMERRAdata

dtypeNetCDF

ydef  144linear-901.25

xdef  288linear-1801.25

zdef  21levels1000975950925900875850825800775750725700650600550500450400350300

tdef  396linear00Z01JAN19791mo

vars3

qv21t,z,y,xSpecifichumidity

u21t,z,y,xEastwardwindcomponent

v21t,z,y,xNorthwardwindcomponent

endvars

然后open***.ctl就行了，之前的问题是打不开ctl文件，怎么改也不行，后来换了台机子就好了。

。

。

所以我说我也不知道怎么回事，希望对你有帮助。

我之前就这样做的，能打开ctl文件，但是d之后，都显示allundefinedvalues，我的ctl如下，麻烦帮我看看哪里错了？

dset  ^%y4%m2.nc

undef1e+15

optionstemplate

titleMERRAdata

dtypenetcdf

xdef  288linear-179.3751.25

ydef  144linear-89.3751.25

zdef  21levels1000975950925900875850825800775750725700650600550500450400350300

tdef  396linear00Z01JAN19791mn

vars3

qv21t,z,y,xSpecifichumidity

u21t,z,y,xEastwardwindcomponent

v21t,z,y,xNorthwardwindcomponent

endvars

可以的，但文件名一定要连续

这个时间长度一定要和你的文件对应好

最近发了一个利用grads批量合并nc文件的帖子  很多人也都问怎么批量描述，，很多人也都去自己尝试，尝试过程中出了各种问题。

推荐兰溪给nc文件写ctl的帖子

“SDFfilehasnodiscernableXcoordinate”的同学们，非常推荐这个帖子。

肯定很多人已经看过这个帖子，也照着做过，还是出了问题。

其实论坛还有很多其他有关ctl批量描述nc的帖子，大家参考一下，对照自己出现的问题，应该大部分情况下都能解决。

但是有一些人很着急，连一个nc文件的描述文件都写不对，就直接批量的，那肯定只会出更多的错误，还一时半会儿改不对。

只会造成时间的浪费。

  我不是来说教的，只是我觉得什么事儿都是循序渐进的，不要那么浮躁。

俗话说磨刀不误砍柴工，其实真的是这么回事。

就说批量描述nc文件，它也是在正确描述一个nc文件的基础上来修改的。

只要描述一个的出来了，往上加一两个语句，用正确的替代格式替代了文件名，再改一下时间长度就出来了。

那能正确描述一个nc文件的ctl就是你磨得很锋利的砍柴刀，有了它后面的柴就相当好砍了。

  下面我就以ncep逐日再分析资料为例说一下nc文件的描述和批量描述。

文件名连续，hgt.1948.nc

hgt.1949.nc  hgt.1950.nc`````````

首先看一下它自带的ctl，如下图红色圈起来的部分，就不详细解释了，对grads有一定了解的人，看一下就知道什么

自带ctl.PNG（29.64KB,下载次数:

9）

下载附件 保存到相册

2013-5-2222:

42上传

下面就可以根据这个自己编写一个ctl文件来描述这个nc文件了。

  照着自带的写下来先，存成1948.ctl然后使用open命令打开，画图，发现出来的图和缺测值设置错误时差不多。

那就想方设法的改缺测值，比如用在grads中使用qattr查看有一句missing_value32766，修改缺测值，再画图，还是那样；qundef出来的是-999000000.000000，再改，再画还是那样。

无论怎么改，都还是那样。

  如果你这么反复折腾，最后还没发现问题，那就说明你没有好好利用论坛的搜索功能，也没有看到兰溪的帖子，也没看到黎大叶子的用Fortran批量为nc写ctl的帖子

  但是不用着急，我看见了，其中最重要的是打开netcdf格式数据的描述文件是需要用xdfopen命令的，那就先要去看看xdfopen能打开的ctl需要怎么写http:

//grads.iges.org/grads/gadoc/gadocindex.html。

  看完了，差不多就能明白了，有这么多前人的成果，那就照着修改呗，最终我修改出来了，图画的也很正常了。

写出来的如下

1.dsetF:

/ncep/daily/hgt.1948.nc

2.titlemeandailyNMCReanalysis（1948）

3.undef-999

4.xdeflon144linear02.5

5.ydeflat73linear-902.5

6.zdeflevel17levels10009258507006005004003002502001501007050302010

7.tdeftime366linear00Z01Jan19481440mn

8.vars1

9.hgt=>hgt17  t,z,y,x  meanDailyGeopotentialheight

10.endvars

复制代码

  光看着正常不行啊，需要和原始的图对比验证了才能确定是对的吧。

所以在grads里面用sdfopen命令打开hgt.1948.nc画第一层第一个时次的图，再用xdfopen打开你编写的ctl，也画第一层第一个时次的图看看。

我擦！

！

！

！

对不上啊，看来还是有问题啊。

说明上面的ctl是有问题的，还得改进才行。

  其实到这步已经基本成功了，仔细看看叠加起来的那两张图，几乎是一样的，只是南北的方向是反的。

那就好办了，在ctl里加上optionsyrev，告诉grads南北要反向不就行了么。

于是最终的ctl出来了，如下

1.dsetF:

/ncep/daily/hgt.1948.nc

2.titlemeandailyNMCReanalysis（1948）

3.optionsyrev

4.*yrev表示y轴反向

5.undef-999

6.xdeflon144linear02.5

7.ydeflat73linear-902.5

8.zdeflevel17levels10009258507006005004003002502001501007050302010

9.tdeftime366linear00Z01Jan19481440mn

10.vars1

11.hgt=>hgt17  t,z,y,x  meanDailyGeopotentialheight

12.endvars

复制代码

再用xdfopen打开画图，这回就一模一样了啊。

说明成功了！

那下面的批量描述就太简单了，比如我要批量描述1948和1949两年的，算一下一个闰年一个平年，一共有时次366+365=731，那么就修改ctl吧

dsetF:

/ncep/daily/hgt.%y4.nc

titlemeandailyNMCReanalysis

optionstemplate

optionsyrev

*yrev表示y轴反向

*undef-999

xdeflon144linear02.5

ydeflat73linear-902.5

zdeflevel17levels10009258507006005004003002502001501007050302010

tdeftime731linear00Z01Jan19481440mn

vars1

hgt=>hgt17  t,z,y,x  meanDailyGeopotentialheight

endvars

就这么简单，只要用%y4代表四位年，把总的时次改成731，再增加一句optionstemplate就可以了。

然后就可以利用xdfopen命令打开画图和原来的对比一下，一样一样的吧，可以批量描述了吧！

注意：

有人可能注意到我把缺测那一项给注释掉了，其实这个是完全不需要的，去掉或者改成任何值，都不影响。

这是问什么？

不细说了，因为我没仔细研究（比较晚了，我要去睡觉了，以后再说，或者有人知道可以告诉我）

其他的资料都同理了，思路都差不多，变通一下就可以了。

1.首先编写可以正确描述一个资料的ctl，并能正确出图。

2.修改该ctl里dest后的文件名，使用合适的替代格式替代（具体的格式看实用手册）

3.增加optionstemplate，这个是批量描述必须加上的一句

4.计算出所要批量描述的文件的总时次，修改ctl里的总时次

5.画图进行验证，如果不对，再根据具体情况做出相应的修改

6.大功告成，可以用了

7.需注意的一点就是optionsyrev，这是可选项，要根据资料实际情况来使用。

今天这个帖子有点儿罗嗦了，大家捡有用的看就行。

其实说那么多，无非是想和一些人说我们这个年代的人已经是站在前人的肩膀上了，很多东西别人都已经有了经验，分享给你你就要好好利用，不要把时间都浪费在发帖和等回复上面。

很多东西，只要自己有一定的基础，在一定合理的范围内就会找到那个答案的·····

[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'

'setgxoutfwrite'

'setfwriteg:

\mls\1.20.32.dat'

t5=1

*作文件名循环

while（t5<=20）

'sdfopeng:

\mls\h2o.'t5'.nc'

      'setx1144'

      'sety191'

      'sett1'

      z=1

      while（z<=25）

      'setz'z

      'dh2o.'t5

      z=z+1

    endwhile

t5=t5+1

endwhile

'disablefwrite'

grads提取某个范围内的资料时注意的一个问题

热度4已有524次阅读2011-11-1916:

10|资料

这两天被一个资料的提取程序弄的晕乎乎的，今天终于发现问题所在了。

提取资料的时候，经度的范围即使原来就是全球的，也不能写成setlon0360，得写成setlon0357.5，就这个2.5（格点间距）把我折腾了这么久，下次一定要注意了！

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

报错：

Datarequestwarning:

requestiscompletelyoutsidefilelimits

gs如下：

'reinit'

'setgxoutfwrite'

'setfwrited:

/sst/4-6nino.grd'

yy=1965

while（yy<=2013）

'sdfopend:

/sst/ersst.'yy'04.nc'

'sdfopend:

/sst/ersst.'yy'05.nc'

'sdfopend:

/sst/ersst.'yy'06.nc'

;'setz1'

;'sett1'

;'setlat-55'

;'setlon170240'

'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'

'ds'

'close3'

'close2'

'close1'

yy=yy+1

endwhile

'disablefwrite'

;

我编写的批量描述文件如下：

DSETd:

/TRMMdata/3b42/9801/3b42.2000%m201.00.7a.nc

TITLETRMMdataofyear2000

DTYPEnetcdf

OPTIONStemplate

UNDEF-9.99E33

XDEF144LINEAR02.5

YDEF20LINEAR02.5

ZDEF1LEVELS1000

TDEF12LINEAR01JAN20001mo

VARS1

pcp=>pcp0t,z,y,xprecipitation

ENDVARS

gs文件如下：

'reinit'

'opend:

/trmmdata/work/2000.ctl'

'setlat060'

'setlon030'

'setlev1000'

'dave（pcp,t=1,t=12）'

'reinit'

 用GrADS转换nc数据

                                    ——byArtmunichfrom

  kittyhare在一篇帖子里向我们初步介绍了使用GrADS转化nc为grd格式（摸我去看看~），我想能不能再详细点，于是借鉴以前自己见到的一个教程，更详细的介绍一下这方面的知识，希望大家指正。

  先给出一个单层的二进制文件转化的gs

'reinit'

'sdfopenf:

/data/nc/1980/air.1980.nc'

'setgxoutfwrite'

'setfwritef:

/data//air.1980.bin'

'setlon0357.5'

'setlat-9090'

'setlev1000'

'sett1640'

'dair'

'reinit'

  但我们知道nc文件是按照经度、纬度、高度、变量、时次顺序排列，要转nc文件，需要靠循环，由于高度的不连续性，我们可以在时间循环里面把每层的高度写出来。

'reinit'

'sdfopenf:

/data/nc/1980/air.1980.nc'

'setgxoutfwrite'

'setfwritef:

/data//air.1980.bin'

'setlon0357.5'

'setlat-9090'

t=1

while（t<=365）

'sett't

'setlev1000'

'dair'

'setlev925'

'dair'

*Youcancontinuetowritehgt

t=t+1

endwhile

'reinit'

  当然，如果用z坐标系，这样z就是从1到17的，在时间的循环里嵌套z=1;while（z<=17）我认为也是可行的，不过自己没试过，有用过的朋友可以告诉一下结果。

高度循环的：

while（t<=8）

'sett't''

z=1

while（z<=21）

'setz'z''

'setlon90125'

'setlat1035'

'definet=tmpprs'

'definerh=rhprs'

'definees=（6.112*exp（17.67*（t-273.15）/（t-29.65）））'

'defineq=rh*（0.62197*es/（lev-es））/100.'

'definee=lev*q/（0.62197+q）+1e-10'

'definetlcl=55.0+2840.0/（3.5*log（t）-log（e）-4.805）'

'definetheta=t*pow（（1000./lev）,（0.2854*（1.0-0.28*q）））'

'definethetse=theta*exp（（（3376./tlcl）-2.54）*q*（1.0+0.81*q））'

'dthetse'

z=z+1

endwhile

t=t+1

endwhile