Fortran+GrADS站点作图详解
mofangbao@气象家园
帖子导读：
以前论坛也有过相关的帖子，但是到目前为止，仍然不断的有朋友看了那些帖子之后还是不知道如何下手，于是我打算为纯粹的新手写一份说明，这份文档会非常详细，因此适合刚接触GrADS的朋友来参考。
帖子介绍了fortran文件中读写数据的基本知识，grads的站点文件的存放格式，如何通过fortran来写出这个格式的文件，如何建立站点文件的CTL文件，如何为站点数据生成映射文件，站点映射文件的作用，如何通过站点CTL文件来预览站点数据，如何用fortran生成grads支持的格点文件，如何给格点配置CTL文件，如何编写GS文件来进行插值，如何绘制等值线以及阴影图，如何屏蔽区域外的图形等内容。如果你对以上问题存在疑惑，希望这个帖子能够对你有所帮助，因此，看完这份文档之后，您应该可以做到这些：
（会的就跳过吧，我也是犹豫了很久要不要写的，算是新手礼包的一部分吧）
开始之前，先了解站点数据的基本内容
一般情况下，我们看到的站点数据是这样的：
50353 126.39 51.43 100 89
50632 121.55 48.46 100 54
50527 119.45 49.13 100 29
50434 121.41 50.29 100 46
50557 125.14 49.10 100 83
50745 123.55 47.23 100 80
50756 126.58 47.26 100 74
50788 131.59 47.14 100 36
50873 130.17 46.49 100 31
50978 130.57 45.17 100 59
第一列为站号，第二、三列为经纬度，第四、五列为其他数据。当然，也可能是把经纬度单独存放，然后数据单独存放，比如气候中心免费开放的160站的月平均温度降水值。如果是micaps资料可以查看micaps用户手册。这个例子中使用的就是micaps的第三类数据格式，示例文件可以从这里下载：
（资料下载）
除了了解这些，还必须了解该数据的缺测值用什么来表示，这个在后面的数据描述文件中要用到。
1、在fortran中读入这些资料
在使用fortran之前，请确保你已经安装了fortran编译器，如果你是xp系统，建议安装CFV6，如果是win7系统，那么可以安装microsoft fortran powerstation4.0（以下简称4.0），这两个软件在气象资料站均有下载，安装完成后，请在你磁盘的某个文件夹下面新建一个文件夹用于本次作图，如果使用的是CVF6，该文件夹的路径中不要包含中文名。新建文件夹完成后，从开始菜单打开fortran编程窗口，然后点击file->new，新建一个free format的自由格式fortran文件,4.0的直接建立一个TextFile即可：
CVF中新建fortran文件
4.0中新建TextFile即可
注意：CVF中先选路径和文件名，然后会建立该文件，4.0中直接建立一个临时文件，当你按保存按钮时才会提示你选择保存的路径，如果你之前用4.0的fortran 编译过其他文件，那么请先关闭工作区再新建，从file->close workspace可以关闭。
当然也可以直接建立一个文本文件，然后重命名为sta2grd.f90，记住，要修改后缀为.f90，表示该文件为自由格式的fortran源文件。现在就可以开始着手读取数据了。
（1）先把fortran的开头和结尾写了，并且用implicit none声明不能使用未定义的变量，同时添加一个注释，自由格式中可以用一个英文的!来表示该行该符号往后均为注释内容，不参与编译：
Program sta2grd
Implicit none
!这里是程序的变量声明
!变量声明结束
!程序开始
!程序结束
End
（2）定义几个变量
Fortran中常用的变量类型有整型：integer，单精度实型real，双精度实型double，字符型character，关于fortran变量的更多内容请参考彭国伦的fortran95程序设计。
通过观察我们要读取的文件：
diamond 3 2011年05月降水量值
0 0 0 0 0 0 1 0 0 1 1000
50353 126.39 51.43 100 89
50632 121.55 48.46 100 54
…
第一行和第二行对我们转换没有什么作用，因此我们只需要定义必须要用到的变量，站号stid，类型为字符型，长度为8，其余的经度lon、纬度lat，一个无意义的变量（列中的100）var和降水数据rain，均为实型，原因在下面的GrADS中的存放方式有说明。
Character*8 stid
Real lon,lat,var, rain
（3）打开数据文件并读入数据
Fortran中使用open语句来实现打开文件，同时指定一个打开的文件号和文件路径，以及其他一些指定参数，下面的语句即可打开这次示例的文件：
Open(1,file=’raindata.txt’,status=’old’)
这里的1是一个文件号，是读写时需要用到的，file后的是文件名，直接写文件名表示在和程序同一个文件夹，status表示当前文件的状态，old表示该文件已经存在，这样就能避免通过write语句误写入数据覆盖了数据文件。
（4）开始读入数据到变量中
在fortran中使用read语句来读入数据，每读取一行文件的指针就会自动向下移动一行（关于文件指针指向这一行你可以想象成该行被选中，要从该行读取数据了）比如read(1,*)var，就表示从打开的1号文件中按默认格式读取数据到var中。我们先跳过两行“没用”的数据：
Read(1,*)
Read(1,*)
这样后面不接任何变量就表示把这两行跳过去了，现在文件指针在第三行的开头了，也就是我们需要的数据了，来读一行：
Read(1,*)stid,lon,lat,var,rain
Fortran 会自动把这五个数值分别给后面的五个变量了，现在文件的指针已经在第四行的开头了，我们需要继续读取数据，当然不能够又写一行，应该继续使用上面这一行来读取，可以通过goto语句实现，goto语句的含义是让程序跳转到指定的标号执行，这里就应该是去刚才的那个语句执行，所以在刚才的Read(1,*)stid,lon,lat,var,rain 前面加上一个标号10，然后在下面写上一句goto 10，如下所示：
10 Read(1,*)stid,lon,lat,var,rain
Goto 10
这样，程序就会不断的往下读取了，但是，文件的行数是有限的，这样做一旦到文件结束，fortran就会报出运行时错误，通常是end-of-file…的错误，这时候可以通过END标签来避免这种情况的出现。End标签能够让read语句读取到文件末尾的时候自动跳转到指定的
标号去执行，再次修改这个读取的语句：
10 Read(1,*,end=100)stid,lon,lat,var,rain
Goto 10
100 contine
这样，当文件读取结束时，就会到100 continue来执行了，continue的意思就是接着按顺序往下执行，本身没有什么含义。
我们现在只是读取数据，没有对数据进行任何输出和处理，读取完毕后应该及时关闭这个文件，使用close方法即可。在100 continue的下面加上close(1)，1同样是open里的那个文件号。
这段程序就完成了这个示例文件的读取，每读取一次，上面几个变量的值就会被刷新一次，如果需要验证读取是否正确，可以在标号为10的read语句和goto 10之前加上一个输出到屏幕的命令，然后用pause命令暂停执行，以便检查：
10 Read(1,*,end=100)stid,lon,lat,var,rain
Print*,stid,lon,lat,var,rain
pause
Goto 10
现在，这段程序应该是这样子的：
Program sta2grd
Implicit none
!这里是程序的变量声明
Character*8 stid
Real lon,lat,var, rain
!变量声明结束
!程序开始
Open(1,file=’raindata.txt’,status=’old’)
read(1,*)
read(1,*)
10 Read(1,*,end=100)stid,lon,lat,var,rain
Print*,stid,lon,lat,var,rain
pause
Goto 10
100 continue
Close(1)
!程序结束
End
这时候就可以编译运行一下了，4.0中点击窗口中左上角的build按钮（shift+f8），弹出的对话框中选择是，然后下面的消息窗口会提示:sta2grd.exe - 0 error(s), 0 warning(s)，表示编译通过了，可以开始执行了，选择build菜单的Execute sta2grd.exe，就可以开始执行了；
CVF中直接点击红色的感叹号左侧有两个朝下箭头的按钮编译，无错后点击红色感
叹号可以执行，每行输出后都会暂停，等按下回车后会继续读取下一行并暂停，经过对比，发现读取成功，然后删除用于测试的print*和pause语句：
4.0中的编译
2、按照GrADS的规定将数据进行输出
接下来所增加的语句都在10号语句和goto 10语句之间，因为这时候几个变量的值是刚刚从文件中读出来，所以要及时把他们写出去。
在输出为GrADS可用格式之前，必须先了解GrADS中的站点数据是如何存放的，只有知道GrADS中规定的存放格式，才能够明白输出时该如何写。
GrADS中的站点数据也是二进制格式，并且按照时间来分组，（本例子中仅有一个时次，多时次的相信你看完后也会写出来了），每个时次下面包括了N个站该时次的数据，每个站的信息由一个表示该站的文件头信息以及数据值组成，每个时次结束后写入一个特殊的标记表示这个时次结束，文件头信息结构如下：
struct reportheader {
char id[8]; /* 站号 */
float lat; /* 纬度*/
float lon; /* 经度 */
float t; /* 相对于该时次的相对单位值 */
int nlev; /* 总层次 */
int flag; /* 表示有无高空资料 */
};
上面的定义是C语言中的定义，转换到fortran中就是我们开始定义的几个
变量，但是还没有定义t(避免和总时次混淆，定义为tim)、nlev和flag，所以先去上面的定义变量的地方补上这上个变量，把real那一行改为：real lon,lat,var, rain,tim ，在real行下面增加一行：integer nlev,flag。这几个变量的详细说明如下：Id:表示站号，每站一个，不能重复，必须以8字节的字符型变量表示（程序中的stid）；
Lat和lon分别表示纬度和经度，为实型（C中的浮点型）这里的顺序和我们示例文件中是反的，因此写入的时候要注意顺序；
T:实型，当有高空资料的时候，使用t来表示一个偏移，比如有整点航空报的时候，如果刚好表示的是整点（比如12点）的数据，那么t=0，如果表示的是12点15分的数据，那么t就是0.25（也就是1/4小时，但可能你的数据描述文件中是12点，这个0.25就是一个偏移），t的值是从-0.5到0.5变化的；如果是地面资料，就让它为0吧。因为没有做过相关工作，因此理解的也不够透彻，麻烦高手指出。
nlev：整型，表示总层次，地面变量算1层，如果有高空变量的话，再加上高空的层次数量，如果nlev为0，表示该时次的结束，也就是上面说的一个特殊的标记；
flag：整型，用来表示有无地面数据，如果有则为1，如果无则为0
根据上面的描述，来对照一下我们的示例文件，stid、lon、lat和数据值已经定义好了，并且能够不断的有数据读出来，因为没有高空数据，所以tim=0.0，nlev=1，有地面数据，所以flag=1，因为这几个值不会发生变化了，所以我们把它写到程序的开头去，免得每次循环读取的时候都要赋值，写在open的上面即可：
!程序开始
tim=0.0
nlev=1
flag=1
Open(1,file='raindata.txt',status='old')
….
在写入数据之前，得新建一个文件，文件号为2，不要和之前的重复，使用fortran的open命令，并且把status参数指定为replace即可，这样如果不存在，则新建，如果存在，则覆盖（如果指定为new的话，存在时就会提示出错了），因为grads中需要的是二进制格式，所以加上format参数，并指定为binary，这句话加在现有的一个open之后即可，如下所示：
Open(1,file='raindata.txt',status='old')
open(2,file='sta.grd',status='replace',form='binary') !新增的新建一个二进制文件语句read(1,*)
因为这个文件只要建立一次即可，所以只能加在标号10的语句之前。
现在就可以在10号语句和goto 10之间加上写出到2号文件的语句了，由于是二进制无格式文件，所以write命令后面无需增加格式控制参数，直接“write(2)变量”的形式即可。按照前面展示的那个文件头+数据值的方式，输出语句当然也就是先写站号，然后是纬度、经度、tim值、nlev值，flag值，最后是一个数据值，也就是下面这句话了：Write(2)stid,lat,lon,tim,nlev,flag,rain
就这样，程序会不停的读一行，写一行，一直到文件读取结束，跳转到100号的那一行继续执行，100号下面就是close(1)了，现在可以关闭1号文件了，但是还不能关闭2号文件，还记得上面说的那个特殊标记么，我们这个时次的文件都写结束了，现在应该写入那个特殊标记告诉GrADS，这个时次结束了。然后才能关闭2号文件，所以在close(1)之后加上：nlev=0
Write(2)stid,lat,lon,tim,nlev,flag
Close(2)
现在你的程序应该是这样子的：
Program sta2grd
Implicit none
!这里是程序的变量声明
Character*8 stid
Real lon,lat,var, rain,tim
integer nlev,flag
!变量声明结束
!程序开始
tim=0.0
nlev=1
flag=1
Open(1,file='raindata.txt',status='old')
open(2,file='sta.grd',status='replace',form='binary')
read(1,*)
read(1,*)
10 Read(1,*,end=100)stid,lon,lat,var,rain
Write(2)stid,lat,lon,tim,nlev,flag,rain
Goto 10
100 continue
Close(1)
nlev=0
Write(2)stid,lat,lon,tim,nlev,flag
close(2)
!程序结束
End
如果看完上面这些，现在你至少应该大概明白GrADS中的站点格式文件是如何存放的了，以及如何写出单时次，单变量的二进制文件。
注意：
A、关于多时次
多时次的需要读取不同时次的变量，然后还是按照这样的格式写入，每个时次写入完成后都应该加一个结束标记，但是多时次的我还是建议你不要第一行写一行，因为那样控制时间的循环对于初学者可能有点复杂，可以先建立一个二维的数组rain(tt,sta)，其中tt是总时次，sta是总站点数，然后经纬度以及站号也相应的改为数组lon(sta)，lat(sta)，stid(sta)，这样等所有数据都读取完毕后再开始写出到二进制文件，就可以这样做：
!这段程序不是本次例子的内容
Do i=1,tt
tim=0.0
nlev=1
flag=1
do j=1,sta
Write(2)stid(j),lat(j),lon(j),tim,nlev,flag,rain(i,j)
enddo
nlev=0
Write(2)stid,lat,lon,tim,nlev,flag
Enddo
Close(2)
这样先读出来的缺点就是你在生成数组之前必须知道你要读取是数据有多少行，这样才能定义数组或者使用动态数组，而且这样做也更耗费内存（虽然现在可用内存都很大了）。你可以试着将上面的程序改成这中写法的，然后在下面贴出来会有分加哦。
B、关于多变量，多层次
如果有多个变量，分两种情况，如果是地面数据，那么直接在rain的后面接下一个变量即可比如
Write(2) stid(j),lat(j),lon(j),tim,nlev,flag,rain(i,j),temp(i,j)
这样就写出了降水值和温度值。
如果是高空变量，要接在地面变量后面写，比如
Do i=1,tt
tim=0.0
nlev=11 !注意，有了高空的，这里的值是10个高空层+1个地面层
flag=1
do j=1,sta
Write(2)stid(j),lat(j),lon(j),tim,nlev,flag,rain(i,j),temp(i,j)
!接下来写高空的10层
do k=1,10
write(2)lev(k),z(i,j,k) !这里的lev和z也都要定义为real型
enddo
enddo
nlev=0
Write(2)stid,lat,lon,tim,nlev,flag
Enddo
Close(2)
C、缺测值的处理
如果资料中包含了缺测值，而且缺测值不是用一个数值来表示的，比如把因为不同原因造成的缺测值用不同的值来表示，那么在write之前要判断是否为缺测，如果是缺测，要统一为同一种缺测值，比如：
do j=1,sta
if(rain(i,j)=-32766.or.rain(i,j)=32766)then
rain(i,j)=-32766
endif
Write(2)stid(j),lat(j),lon(j),tim,nlev,flag,rain(i,j)
Enddo
到现在，已经将数据写成了grads可以识别的二进制格式了，如果没有修改输出文件路径，sta.grd就在这个f90程序的当前路径下。
3、为生成的站点数据编写CTL（数据描述文件）
由于这个生成的数据并不是很完整的自描述文件（比如可以用sdfopen命令打开的一些nc文件），所以，还得给这个数据配备相应的数据描述文件，这样GrADS才能知道这是一个站点格式文件，以及这个文件中有多少时次，有没有高空层次变量，多少变量等信息。总的来说，关于数据描述文件的更多内容，可以参考grads的官方手册，或者lasg编写的中文手册。这里仅仅介绍和此次例子有关的站点格式的数据描述文件，一个标准的站点格式数据描述文件应该如下所示：
DSET D:/station_data_sample.dat
DTYPE station
STNMAP D:/station_data_sample.map
UNDEF -999.0
TITLE Station Data Sample
TDEF 10 linear 12z18jan1992 12hr
VARS 11
ps 0 99 Surface Pressure
ts 0 99 Surface Temperature
ds 0 99 Surface Dewpoint Temperature
us 0 99 Surface U-Wind
vs 0 99 Surface V-Wind
elev 0 99 Elevation of Station
z 1 99 Height
t 1 99 Temperature
d 1 99 Dewpoint Temperature
u 1 99 U-Wind
v 1 99 V-Wind
ENDVARS
DSET用来说明二进制的站点数据的存放位置；
Dtype station用来说明这个二进制数据是一个站点格式的数据，要按照站点格式的数据规范来读取；
STNMAP 后面的文件是一个站点映射文件，这里指定的是一个含路径的文件名，具体的等介绍完CTL文件后再介绍，它能够加速GrADS对站点数据的处理；
UNDEF后面接的是缺测值，也就是数据中所包含的缺测值，这个在第一步了解数据格式的时候就应该明确自己的数据缺测值是什么，并且在fortran输出的时候做统一的处理。
TITLE是这个数据的文件名，可以自己指定
TDEF是这个数据有多少时次，10表示有10个时次，linear表示时次的增长是固定间隔的，后面接的时间是第一个时次的时间，上面的ctl中是1992年1月18日12点，12hr表示下一个时次是12小时之后。
VARS 11表示有11个变量，也就是多变量的，这里要注意的是下面的变量顺序必须和你用fortran写出的时候的变量顺序完全一致，先地面层变量，后高空层变量。Ps 0 99 surface pressure就表示第一个写入的变量定名为ps，0表示它没有高空方向上的数据，99很多人不理解，对于一般的二进制文件，输入99即可，因为这个值是被忽略的，但是必须要写，但是对于dtype为其他格式的文件时，就需要按情况填写了，具体的可以参考grads的官网，（点我前往介绍页面）。后面的z 1 99 Height表示z为一层高空变量。
ENDVARS表示变量说明结束。
根据上面的介绍，就可以写出这个示例文件的CTL了，需要改变的除了文件路径，主要的就是时次只有1层，变量没有高空层，只有一个地面层的，修改如下：
DSET E:/program/example/sta2grd/sta.grd
DTYPE station
STNMAP E:/program/example/sta2grd/sta.map
UNDEF -32766
TITLE Rain Data Sample
*这里的时间可以根据上面示例文件中的第一行的时间填写
TDEF 1 linear 01may2011 1mo
VARS 1
rain 0 99 Rain Data
ENDVARS
现在把这个ctl文件保存在和sta.grd相同的目录吧，我的是
E:/program/example/sta2grd/sta.ctl （*和#在ctl和gs文件中可以注释一行，且必须放在一行的开头，前面不能有空格等其他任何字符）
4、为站点的二进制数据生成站点映射文件（*.map文件）
站点映射文件的具体细节也许不太好理解，你只需要知道他是为了让GrADS 能够能够快速的处理大量站点数据而生成的，生成这个文件必须运行一个外部的命令叫做 STNMAP 命令，这个文件和grads.exe是在同一个路径下的，如果你对CMD命令中如何切换目录到grads.exe那个目录不太理解，或者你不想去输入grads.exe所在的路径，那么可以直接在grads打开的情况下运行下面的命令：！stnmap –i E:/program/example/sta2grd/sta.ctl
正常情况下会出现类似下面（不同版本有所不同）的一些信息：
请注意，在使用stnmap 命令的时候前面有一个英文的感叹号，这表示在grads窗口中运行外部命令，因为grads的当前路径和stnmap是在同一个路径的，所以不需要在stnmap 前面加上stnmap.exe所在的路径了；stnmap后面的路径是sta.ctl的路径，路径中请用/，
而不是\。
这时候，在sta.ctl所在的目录就会生成一个叫做sta.map的文件，这个文件的文件名就是在sta.ctl中stnmap行指定的那个文件（在ctl指定的时候该文件还没有生成呢，如果你指定到其他文件夹去了，现在也得去其他文件夹找了）。
提醒：生成的.map文件和当前的路径有关，如果你生成之后，改变了有关的文件路径，需要重新生成.map文件，如果是拷贝到另一个电脑上，也建议重新生成（不过测试发现如果两台机器的所有文件路径及文件名完全相同，.map文件可以重用）。
5、看看站点数据如何
现在有了数据，有了ctl文件和map文件，就可以通过grads查看站点文件的数据。
打开grads窗口，输入open e:\program\example\sta2grd\sta.ctl 将这里的路径换为你的路径，然后grads会显示出相关的维数信息，接着进行简单的经纬度设置，set lon 70 140 ,set lat 15 55，让一会儿的数据只显示中国范围的底图，然后d rain这个ctl中的变量名，像下面的窗口所展示的那样：
最后通过一句printim命令输出图形，后面的white参数表示用白色作为背景色，输出的图像如下：
现在grads已经把示例的站点数据按照经纬度输出了。
6、生成一个格点的背景场，为插值做准备
因为需要输出的是等值线图、阴影图等用格点文件才能够画出来的图，所以，需要把站点文件再插值到格点上，这样就能借助grads强大的绘图命令绘制各种图形了。
小知识：所谓的格点文件，就是像小时候用过的田字格一样，在每个行列交叉的地方都是一个格点，在grads中也就是经纬网交叉的地方，现在的站点文件要插值为格点，就必须要有一个现成的经纬网来做参考，让grads找到那些交叉点进行插值，至于原来交叉点上的值不用知道，用到的只是它交叉点的坐标信息而已，关于插值的细节，grads会自动处理，不需要我们知道。因此，现在要做的就是去生成一个经纬网，在这里称为格点背景文件。
在制作该文件之前，需要明白自己要插值的分辨率是多少，这样才能在fortran程序中对分辨率进行控制。
分辨率简单的说就是间隔，比如纬向分辨率，就是东西方向上360度分为若干个个点之后，每个格点之间的间隔。如果分辨率为1*1那么全球就有360*181个格点。如果你需要2.5*2.5的分辨率，全球就有144*73个格点。
现在假设需要的是2.5*2.5的中国区域的格点，从东经70到140度，北纬15到55度，
这样就可以划分出29*17，共493个格点。那现在就可以去写fortran程序来生成格点数据
了。
写到这里才忽然想起来，原来写站点画图的参考手册也必须要涉及格点数据的生成方法，真是一举两得了。
（1）新建fortran90的源文件
按照刚开始时候介绍的步骤，打开fortran后，先关闭之前打开的工作区（4.0的才要这么做，CFV6重新打开时默认是没有工作区会打开的），然后建立一个名为grid.f90的fortran 文件，路径还是在sta2grd这个文件夹下面。然后把程序的框架先写出来：Program gen_grid
Implicit none
!这里是定义的变量
!程序开始
!程序结束
End
（2）定义变量，写出格点场
前面说过，这次需要的是中国区域的29*17个格点，所以定义两个整型常量x、y分别表示经纬向的格点数：
Integer,parameter::x=29,y=17
因为implicit none定义了不能使用未定义的变量，所以要为下面的循环中使用到的i和j这两个整型变量进行定义：
Integer i,j
因为数据文件是一定要有数据的，所以定义一个变量grid,用来表示那个数据，因为前面说过，在插值中，这个数据是用不到的，只用到它提供的格点场，所以就让数据值永远为1好了，还记得前面说过grads中的数据都是浮点型，也就是fortran中的单精度实型，所以可以这样定义：
Real grid
然后接着将其赋值为1
grid=1.0
接下来打开一个二进制的新文件grid.grd，用来存储值为1的格点场：
Open(1,file=’grid.grd’,status=’replace’,form=’binary’)
现在这个1号文件grid.grd已经打开，可以按照grads格点数据的存放格式往里面存放数据了。
GrADS中格点数据是按照x y z var t的顺序来存放的，也就是说第一个变量是第一时次，第一层次，第一个纬度上，第一个经度格点的值，第二个变量就是其他不变，第二个经度格点的值(x+1了)，等所有经度的格点（比如360个）放完了，纬度再挪动一个(y+1)，然后再所有经度的格点，等纬度挪动结束（比如全球的91个），高度层向上挪一层(z+1)，等高度层结束了，再把时次变为第二个时次(t+1)，所以，在多个变量存放的时候应该是时间的变化在最外层，然后存放x y z上的第一个变量，而x y z 也是按照z循环在外，内层为y 循环，最内层才是x循环的顺序。
因此，代码应该如下考虑：
Do i=1,y
Do j=1,x
Write(1)grid
Enddo
Enddo
写完后不要忘了关闭文件：
Close(1)
说明一下，29*17是如何算出来的，因为间隔是2.5，所以(140-70)/2.5后取整再+1即可（注意0到360度有一个点是重合的，所以不用加1了，就是144个格点）。
在写格点场用来插值的时候，只要很简单的这样写就行了，不需要考虑多时次，多层次的问题，在插值的时候永远指定使用t=1时候的值就行了，即便你的站点数据是多时次的。
顺便指出，将文本格式的格点数据写出为grads支持的二进制格式时用来作图（而不是仅给站点插值）的时候，需要考虑t和z，要按照前面说的顺序来进行循环，比如下面的代码模板：
Do i=1,tt !最外层的时次
Do j=1,zz !第一个变量的层次循环
Do k=1,y !第一个变量的纬度循环（南北向）
Do m=1,x !第一个变量的经度循环（东西向）
Write(1)var1(I,j,k,m)
Enddo
Enddo
Enddo
Do j=1,zz !第二个变量的类似循环
Do k=1,y
Do m=1,x
Write(1)var2(I,j,k,m)
Enddo
Enddo
Enddo
Enddo
现在你的fortran程序应该是下面这样的：
Program gen_grid
Implicit none
!这里是定义的变量
integer i,j
integer,parameter::x=29,y=17
real grid
!程序开始
grid=1.0
Open(1,file='grid.grd',status='replace',form='binary')
Do i=1,y
Do j=1,x
Write(1)grid
Enddo
Enddo
close(1)
!程序结束
End
编译运行吧。
看一下sta2grd的文件夹下面已经生成了grid.grd这个文件了。
7、为生成的格点场配置CTL文件
前面配置过了站点数据的CTL文件，格点的CTL文件要比站点的多一些内容，格点CTL 文件大概包括以下几个内容：
DSET gridded_data_sample.dat
TITLE Gridded Data Sample
UNDEF -9.99E33
XDEF 180 LINEAR 0.0 2.0
YDEF 90 LINEAR -90 2.0
ZDEF 10 LEVELS 1000 850 700 500 400 300 250 200 150 100
TDEF 4 LINEAR 0Z10apr1991 12hr
VARS 4
slp 0 99 sea level pressure
hgt 10 99 heights
temp 10 99 temperature
shum 6 99 specific humidity
ENDVARS
可以看出，DSET（格点文件路径）、TITLE（标题）、UNDEF（缺测值）、TDEF（时间维定义）、以及VARS以下的变量定义和站点文件的是一致的，变量slp是只有地面层（Z=1时有值），hgt和temp是Z方向有10层的变量，shum在Z方向只有6层（Z=1,6，也就是1000到300）。
不同的地方是多出了XDEF、YDEF、ZDEF这三个定义，这三个是格点数据描述文件所必须的，分别定义x、y、z方向上的纬度变化情况，比如上述的这个CTL文件，在x方向上有180个格点，从0.0开始，格点之间的间隔为2。Y方向也是类似，共90个格点，从-90开始，间隔为2。Z方向没有使用linear来说是按间隔增加，而是使用LEVES直接列出了所有的层次分别是Z=1时为1000，一直到Z=10的时候为100，不等间隔。
按照用fortran写出时的格式，就可以配置这个CTL文件了，先在sta2grd目录下面新建一个名为grid.ctl的CTL文件，然后写入以下内容：
DSET E:\program\example\sta2grd\grid.grd
TITLE Grid Data Sample
UNDEF -9.99E33
XDEF 29 LINEAR 70.0 2.5
YDEF 17 LINEAR 15.0 2.5
ZDEF 1 LEVELS 1000
TDEF 1 LINEAR 01may2011 1mo
VARS 1
g 0 99 Grid Data
ENDVARS
8、编写gs脚本进行插值后作出图形
终于快到要画图的时候了，下面还在sta2grd这个目录新建一个叫sta.gs的grads脚本文件，在grads脚本文件中可以对grads变量进行各种操作的，具体脚本文件的介绍可以参考lasg的中文手册。这里只对用到的命令做一些说明。
Grads的脚本命令中，除了do，while，if等脚本关键字可以直接写出来，其他的均需要将设置放在两个上引号（单引号）之间。Grads的脚本中不能出现tab符号，也就是制表符，否则会报错。
打开这个新建立的脚本文件，在首行输入以下内容：
‘reinit’
这个命令表示关闭所有打开的文件，清空图形输出，并还原所有设置为默认设置。
然后依次打开两个ctl文件，grads中是通过打开ctl文件来打开数据文件的，而不是直接打开grd文件：
'open e:\program\example\sta2grd\grid.ctl'
'open e:\program\example\sta2grd\sta.ctl'
这样，格点背景文件就被打开为1号文件，站点文件就被打开为2号文件，在显示变量的时候，如果变量不是1号文件中，要在变量后面加上一个点号，如这里的sta里面的rain 变量在脚本中要写为：rani.2。
接下来先设置一下一会儿要输出的范围，因为是在中国地区，所以设置经纬度到中国的范围：
‘set lon 70 140’
‘set lat 15 55’
Grads中的set lon和set lat命令是用来设置经纬度的，后面的两个参数分别是起始值和结束值，都是按照经纬度来算的值。
为了更清楚的显示中国区域，把地图设置为高精度的地图：
‘set mpdset hires’
其中的set mpdset命令是设置地图的命令，地图存放的地点在不同的grads版本中稍有不同，你可以在grads的安装目录中搜索hires这个文件，找到后，如果以后要新增地图文件，直接放在相同的文件夹就行了（在气象家园打包的grads中已经预置了很多常用的地图了）。
接下来就插值并绘图吧，关键的一句话了：
‘d oacres(g,rain.2)’
这句话会先使用oacres函数进行插值，然后把插值的格点场进行绘图输出，d就是display 的缩写，是grads中的输出图像命令，因为插值后是二维格点数据，所以默认的图形是等值线图。
现在还想把图像输出，继续使用前面介绍过的printim命令：
‘printim e:\program\example\sta2grd\sta_grid.png’
在这句话的下一行加上一个分号;，表示这个gs脚本的结束
现在你的gs文件应该是这样的：
'reinit'
'open e:\program\example\sta2grd\grid.ctl'
'open e:\program\example\sta2grd\sta.ctl'
'set lon 70 140'
'set lat 15 55'
'set mpdset hires'
'd oacres(g,rain.2)'
'printim e:\program\example\sta2grd\sta_grid.png white'
;
现在打开grads软件吧，在软件中输入：
Run e:\program\example\sta2grd\sta.gs
然后回车，图形输出窗口上应该会很快就有图形输出了，同时在sta2grd目录下面也生成了输出的图形：
终于绘出图像了！
9、加上阴影图，以及屏蔽中国区域外的图像
阴影图在grads中通过set gxout shaded来进行设置，但是要注意，阴影图必须在等值线图之前绘制，否则阴影图将会把等值线图覆盖。这次因为要d两次变量，所以可以先用一个变量把插值出来的数据保存，grads中定义变量的命令为define，这个命令只能用来接收格点场数据，而不能是一个数值，在grads脚本语言中，一个数值变量不需要声明即可使用，这些可以参考lasg的中文手册。设置完输出类型后直接d保存的这个变量即可，将上面的d oacres(g,rain.2)修改为：
‘define rgird=oacres(g,rain.2)’
然后进行输出类型的设置，并绘制变量，紧接着这句话下一行增加：
‘set gxout shaded’
‘d rgrid’
‘set gxout contour’
‘d rgrid’
如果你现在用run命令运行你的gs，图形应该是像下面这样的：
是不是感觉有点变扭呢，正如这一小节的标题所说，还想屏蔽中国区域以外的图形，那就要用到一个叫做cnbasemap.gs的文件以及相关的地图文件了，这个我们气象家园论坛的传说超版有专门的帖子说明（点我去看），如果你使用我们论坛打包的grads，就没有这些麻烦了，所有需要的东西都已经打包好并且设置好路径了（点我去下载）。把第二个d rgid（也就是set gxout contour下面的那一行）改成如下语句：
‘cnbasemap rgrid’
为什么要把原来的d去掉，那样岂不是不会输出等值线了？打开cnbasemap.gs文件就会发现在文件里面已经进行了输出了，所以没必要先输出一次，所以现在整个gs文件是这样的：
'reinit'
'open e:\program\example\sta2grd\grid.ctl'
'open e:\program\example\sta2grd\sta.ctl'
'set lon 70 140'
'set lat 15 55'
'set mpdset hires'
'define rgrid=oacres(g,rain.2)'
'set gxout shaded'
'd rgrid'
'set gxout contour'
'cnbasemap rgrid'
'printim e:\program\example\sta2grd\sta_grid.png white'
;
（如果你电脑里面还没有下载cnbasemap.gs建议去下载我们打包的grads 吧，或者自己按照传说的帖子进行配置也可，否则将会不显示阴影图）画出来的图形如下：
如果不喜欢图像下面的grads标识，可以在set lat 的下一行加一句’set grads off’，如果想给这幅图像添加一个色标，那么可以在cnbasemap的下一行加上’cbarn’，加上这两句话之后的脚本如下：
'reinit'
'open e:\program\example\sta2grd\grid.ctl'
'open e:\program\example\sta2grd\sta.ctl'
'set lon 70 140'
'set lat 15 55'
'set grads off'
'set mpdset hires'
'define rgrid=oacres(g,rain.2)'
'set gxout shaded'
'd rgrid'
'set gxout contour'
'cnbasemap rgrid'
'cbarn'
'printim e:\program\example\sta2grd\sta_grid.png white'
;
再次运行，这次出来的图形就好多了：
更多的图形美化设置可以在我们grads版块找到帖子学习。
10、总结
在这个例子中，开始之前先了解了资料的基本信息，开始之后，首先使用fortran输出了二进制的站点数据(sta.grd)，然后为这个站点数据配置了ctl(sta.ctl)，并且还生成了一个map文件(sta.map)链接到这个数据，期间还查看了站点数据的内容，接着根据插值的需要，又使用fortran生成了一个格点场(grid.grd)，之后为这个格点场配置了ctl文件（grid.ctl），总共5个文件，当然另外还有两个f90的fortran源文件以及编译的时候生成的一些文件。在这些工作完成之后，开始编写绘图脚本文件（sta.gs），最终通过这个脚本文件进行了插值，画出了图形，此外还调用cnbasemap文件进行了中国区域外的图形屏蔽。
终于写完了，以上只是我个人经验的一些总结，旨在为接触grads和fortran不久的朋友提供一些所谓的参考，如果你能跟着上面说的步骤一步看下来，做下来，相信用站点画图对你来说只是小事一桩了（当然meteoinfo的界面操作更为快捷^v^）。
如果这个参考对你有了帮助，在你以后grads运用的很熟练的时候不要忘了来家园给新手提供一些力所能及的帮助。
清风2011年12月5日
Grads处理GRIB格式文件的准备 LYanbing 2007-7-2初稿，2007-7-12修改1 开场说明 WHAT IS GRIB? GRIB (GRIdded Binary) is an international, public, binary format for the efficient storage of meteorological/oceanographic variables. Typically, GRIB data consists of a sequence of 2-D (lon,lat) chunks of a (in most general sense) 4-D variable (e.g., u comp on the wind = f(lon,lat,level,time)). The sequence is commonly organized in files containing all variables at a particular time (i.e., 3-D (lon,lat,level) volume). 大气所的NCEP再分析资料使用这种格式。 这里针对6小时一次的1°×1°，26层数据来处理。 Grads中识别路径的方式基本为Unix的方式，即路径中用斜杠/，而不是反斜杠\，cmd中也支持这种方式，所以，使用Grads及其相关组件时，指定路径用斜杠/会很方便。 Grads中用!pwd可以看到当前目录，C:盘对应/cygdrive/c/，d:盘对应/cygdrive/d/。cmd中用pwd看到的也是如此，cygdrive是怎么来的？与cygwin程序有关，它能把Unix程序嫁接到windows下使用，它的目录系统以/cygdrive/为根目录。 如果ctl文件中，数据文件指定不是全路径，而是^，则可以在open命令中指定上述形式的全路径，例如： ga-> open /cygdrive/d/data/ncep/grib2006060100.ctl 它等效于： ga-> open d:/data/ncep/grib2006060100.ctl 2 生成描述文件 PCGrads软件的User’s Guide中介绍了GRIB及其处理方法，但不完全。实际上有两种方法：1）利用工具grib2ctl.exe生成整个文件的描述文件.ctl，再利用工具gribmap.exe生成映射文件.idx；2）利用工具wgrib.exe解码文件中需要使用的部分记录，建立新的数据文件，然后人工建立描述文件.ctl。 为了使用方便，环境变量Path中增加Grads可执行文件所在目录，则在cmd中其他路径下亦可访问所有该目录下的工具。 2.1 方法1 使用工具grib2ctl.exe生成GRIB数据文件的描述文件ctl，之后还要使用gribmap工具生成映射文件.idx。 1. 工具grib2ctl.exe的获得。
Components of a GrADS Data Descriptor File DSET data_filename back to top This entry specifies the filename of the data file being described. If the data and the descriptor file are not in the same directory, then data_filename must include a full path. If a ^ character is placed in front of data_filename, then data_filename is assumed to be relative to the path of the descriptor file. If you are using the ^ character in the DSET entry, then the descriptor file and the data file may be moved to a new directory without changing any entries in the data descriptor file, provided their relative paths remain the same. For example: If the data descriptor file is: /data/wx/grads/sa.ctl and the binary data file is: /data/wx/grads/sa.dat then the data file name in the data descriptor file can be: DSET ^sa.dat instead of: DSET /data/wx/grads/sa.dat If data_filename does not include a full path or a ^, then GrADS will only look for data files in the directory where you are running GrADS. GrADS allows you use a single DSET entry to aggregate multiple data files and handle them as if they were one individual file. The individual data files must be identical in all dimensions except time, and the time range of each individual file must be indicated it its filename. To accomplish this, the DSET entry has a substitution template instead of a filename. See the section on Using Templates for a description of all the possible components of the template. Second, the OPTIONS entry must contain the template keyword. CHSUB t1 t2 string back to top (GrADS version 1.9b4) This entry is used with a new option for templating data files that allows for any user-specified string substitution, instead of only date string substitution. This is useful when none of the standard template options match the time ranges in the files you wish to aggregate, or if the files are located on different disks. When you put the %ch template in your DSET entry, then you also need to put additional CHSUB entries in the descriptor file. The string will be substituted for%ch in the data file name for the time steps beginning with t1 and ending with t2.See the section on Using Templates for examples. DTYPE keyword back to top The DTYPE entry specifies the type of data being described. There are four options: grib, hdfsds, netcdf, or station. If the data type is none of these, then the DTYPE entry is omitted completely from the descriptor file and GrADS will assume the data type is gridded binary.
Fortran+GrADS站点作图详解 mofangbao@气象家园 帖子导读： 以前论坛也有过相关的帖子，但是到目前为止，仍然不断的有朋友看了那些帖子之后还是不知道如何下手，于是我打算为纯粹的新手写一份说明，这份文档会非常详细，因此适合刚接触GrADS的朋友来参考。 帖子介绍了fortran文件中读写数据的基本知识，grads的站点文件的存放格式，如何通过fortran来写出这个格式的文件，如何建立站点文件的CTL文件，如何为站点数据生成映射文件，站点映射文件的作用，如何通过站点CTL文件来预览站点数据，如何用fortran生成grads支持的格点文件，如何给格点配置CTL文件，如何编写GS文件来进行插值，如何绘制等值线以及阴影图，如何屏蔽区域外的图形等内容。如果你对以上问题存在疑惑，希望这个帖子能够对你有所帮助，因此，看完这份文档之后，您应该可以做到这些： （会的就跳过吧，我也是犹豫了很久要不要写的，算是新手礼包的一部分吧） 开始之前，先了解站点数据的基本内容 一般情况下，我们看到的站点数据是这样的： 50353 126.39 51.43 100 89 50632 121.55 48.46 100 54 50527 119.45 49.13 100 29 50434 121.41 50.29 100 46 50557 125.14 49.10 100 83 50745 123.55 47.23 100 80 50756 126.58 47.26 100 74 50788 131.59 47.14 100 36 50873 130.17 46.49 100 31 50978 130.57 45.17 100 59 第一列为站号，第二、三列为经纬度，第四、五列为其他数据。当然，也可能是把经纬度单独存放，然后数据单独存放，比如气候中心免费开放的160站的月平均温度降水值。如果是micaps资料可以查看micaps用户手册。这个例子中使用的就是micaps的第三类数据格式，示例文件可以从这里下载： （资料下载） 除了了解这些，还必须了解该数据的缺测值用什么来表示，这个在后面的数据描述文件中要用到。 1、在fortran中读入这些资料 在使用fortran之前，请确保你已经安装了fortran编译器，如果你是xp系统，建议安装CFV6，如果是win7系统，那么可以安装microsoft fortran powerstation4.0（以下简称4.0），这两个软件在气象资料站均有下载，安装完成后，请在你磁盘的某个文件夹下面新建一个文件夹用于本次作图，如果使用的是CVF6，该文件夹的路径中不要包含中文名。新建文件夹完成后，从开始菜单打开fortran编程窗口，然后点击file->new，新建一个free format的自由格式fortran文件,4.0的直接建立一个TextFile即可：
FNL 1X1 资料认识与应用 兰溪之水整合版2011-12-16 14:34:05（第一版） 2012-05-01 22:42:49（第二版）FNL 1.0X1.0数据下载地址：https://www.doczj.com/doc/9c12703299.html,/datasets/ds083.2/需要注册通过方可下载。 FNL 1.0X1.0 资料包含的物理量及其缩写 缩写参数名称 No4LFTXsfc 近地表四层等压面的抬升指数 No5WAVAprs 500 hPa等压面位势高度距平 No5WAVHprs 500 hPa等压面位势高度 ABSVprs 绝对涡度 CAPE 对流有效位能 CIN 对流抑制能 CLWMRprs 云水 CWATclm 气柱云水 GPAprs 位势高度距平 HGT 位势高度 HPBLsfc 地表行星边界层高度 ICECsfc 海冰密集度 LANDsfc 陆地覆盖 LFTXsfc 地表抬升指数 O3MRprs 臭氧层混合比 POTsig995 位温 PRE 气压 PWATclm 可降水量 RH 相对湿度 SOILW 土壤体积含水量 SPFH 比湿 TCDCcvl 对流云总云量 TM 温度 TOZNEclm 臭氧含量 UGRD u分量 VGRD v 分量 VVEL 垂直速度 VWSH 垂直风切变 WEASDsfc 累积雪量
GrADS处理FNL 1.0X1.0（grib1）数据 处理前需要先将grib2ctl.exe放到GrADS安装文件夹GrADS19\win32（1.9版本）或OpenGrADS\Contents\Cygwin\Versions\2.0.a9.oga.1\i686（2.0版本）下，方便操作。gribmap.exe（这个GrADS自带的） 第一步：先生成一个描述文件ctl 打开命令提示符， C:\Documents and Settings\Administrator>grib2ctl grib_file > grib_file.ctl 如： 或者进入GrADS： gs>!grib2ctl grib_file > grib_file.ctl （！表示调用外部的程序，注意路径用“/”） 如： 这样fnl_20101013_00_00_c.ctl描述文件就生成了！ 第二步：利用GrADS自带的gribmap.exe生成索引文件： C:\Documents and Settings\Administrator>gribmap –v –i grib_file.ctl > b.txt 如： 或进入GrADS： 这样就会在H:\test路径下生成名为fnl_20101013_00_00_c.idx的索引文件，并会把整个映射过程写入到b.txt文件中，当然这里也可省略查看映射过程，即 C:\Documents and Settings\Administrator>gribmap –i grib_file.ctl 这样就可以开始用GrADS读取FNL文件画图了！ 不过我们处理FNL资料的时候一般都是不止一个文件，所以我们就需要用到批处理了！ 第一步：利用grib2ctl.exe生成初始时刻ncep数据的ctl文件；示例中生成的是fnl_20101013_00_00_c数据的文件。
收稿日期:2002-11-1 作者简介:高文良,男,31岁,主要从事短期气候预测等研究工作。 气象绘图软件GrADS 高级使用技巧 高文良1 刘晓燕2 曾小东3 (11成都高原气象研究所 成都 610072; 21阿坝州金川县气象局 624100 31阿坝州马尔康县气象局 624000 ) 摘 要:本文通过介绍气象绘图软件GrADS 使用中 的一些高级技巧,分析了GrADS Script 语言中的难理解之处,并介绍了在GrADS 中做合成分析和t 检验的程序,可以对学习和使用GrADS 的科研人员起参考、帮助作用,达到事半功倍的效果,促进GrADS 软件的应用更广泛、深入,提高科研人员的工作效率。关键词:GrADS ;技巧;程序中图分类号:TP391文献标识码:C 文章编号:1003-7187(2002)04-0057-03 1 引言 气象绘图软件GrADS (Grid Analysis and Dis 2play System )是免费共享软件,可随时从互联网上下载(http ://https://www.doczj.com/doc/9c12703299.html,/grads/),后续版本正陆续推出[1]。GrADS 有丰富的内部函数,可以对数据进行计算和分析处理。它支持处理格点资料和站点资料,并且支持对GriB 码文件、特殊格式文件(如一字节整型、二字节整型、大中型机器二进制数据等)的直接读取,气象科研领域应用非常广泛[2]。在其最新1.8SL 9版本中,GrADS 又将应用领域推进到了海洋学科,功能也得到了进一步地增强和扩展。 但在使用GrADS 的过程中,特别是使用GrADS Script 语言编程当中,有一些问题需要特别注意。因为GrADS Script 语言是一种类似于VB Script (或MA TL AB Script )的高级语言,稍不注意就容易出现错误,且查错十分困难,这点与其他语言有较大的差别。2 站点数据处理 GrADS 中站点数据处理基于格点数据的基础之上,需先将站点数据通过Cressman 客观分析方法内插至格点上,然后再依照格点资料的处理方式对站点资料数据进行分析和处理。所以,内插的背景网格点的选取就显得比较重要。提供背景格点资料数据只起一个背景网格场的作用,格点数据并不参与运算,只提供网格背景,告诉站点数据插至什么点及各个点的距离、综合考虑几个点来插值等信息。 因此,对于需要揭示小尺度特征气象场的站点数据,就需要将背景网格点的间距取得小一些,这样可以将局地小尺度信息完整地体现出来。反之,如果要强调大尺度的信息,就可以将背景格点场的间距取得大一些,将小尺度的噪声滤掉,体现大尺度场的特征。两种方法的最终目的是画出真实而美观的原始数据场的图形。 要将已有的站点数据资料转换成GrADS 可以读取格式的站点数据,也可使用Visual Fortran 6.0或C 语言(Turbo C 、Visual C ++等)进行转换。但必须按照GrADS 的规定格式来变换。数据文件头的结构和定义也必须按规定且与后续数据一致。在Visual Fortran 6.0语言中可以用流式文件(Stream )的方式生成GrADS 需要的站点数据文件,可以写成多时次、多层次和多变量的数据集(在Power Station 4.0中则不能写成多时次和多层次的数据集,它不支持Stream 方式的文件,只能写一个时次的数据)。而在C 语言中的操作更加方便,因为C 语言中对文件的操作方式一般是以字节为单位进行的,没有其他多余的信息。只要按照GrADS 要求的数据格式生成数据,以上的各种语言环境下都能够被GrADS 正确读取。 另一个需注意的问题是在格点文件和站点文件之间,一定使两者的数据描述文件保持时间上的一致性。起始时间和时间步长都要符合数据本身的规定。如果起始时间不一致,在GrADS 中会出现少于两个站的提示,画出的图被标示为缺测值。如果上述两种数据的时间步长不一致,则可能出现错误的结果,导致在一个特定时间下的数据成为另一个时间下的数据,而不是操作者想要的时次。或者图形的结果和前面第一种错误一样,出现全为缺测值的错误情况。
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就可以打开所有文件。 *************************************************************************************************************** *******************
第一讲 GrADS简介 一、G r A D S的应用领域及其功能 *GrADS的全称：“The Grid Analysis and Display System” *应用领域：可在UNIX工作站以及个人微机上进行地球科学领域的数据资料分析和绘图 *功能:对数据进行访问、分析和绘图 1.可以根据需要绘制单线图、直方图、等值线图、填色等值线图、流线图、矢量图、站点模型图等各类图形. 2.用描述语言编程,达到理想的绘图效果 3.可以把在GrADS中绘制的图形以文件的形式保存起来,以备对其进行编辑 4.调用GrADS的内部函数，能够对数据进行某些特定的计算，然后输出计算结果 二、GrADS中常用的基本概念 *数据格式：GrADS能够识别的数据为二进制无格式直接或顺序记录格式，该种格式数据的生成可以通过Fortran语言编程来实现。 *数据类型：格点数据（NCEP/NCAR的再分析资料）、站点数据（站点实测资料）、Grib数据（NMC产品）。*维数环境：GrADS的操作对象为4维的数据集（4D data set）,包括空间三维（纬度、经度、高度）和时间一维。可以固定其中的一维或者几维以获得低于四维的数据子集。此概念是对于格点资料而言的。维数环境的定义可以在两种坐标上进行。一种是地球坐标（world coordinate），以经纬度为度量单位；一种是格点坐标（grid coordinate），以网格点数为度量单位。 *几种文件类型： *.dat―数据文件 *.ctl—原始数据描述文件 *.gs―GrADS控制文件，用命令run执行之 *.exe―GrADS在DOS环境下的各种执行文件 三、启动和退出GrADS *启动GrADS的两种方式： 1.在dos环境下直接输入grads命令 即：切换到MS-DOS方式，进入到 F:\pcgrads\msdos子目录下 输入：grads 回车 2.从windows桌面上的“我的电脑”进入F:\pcgrads\子目录，然后双击g.exe图标 *注意：在启动GrADS时，系统会问你以何种方式进入。此时有四个参数可供选择：b—以批处理的形式运行GrADS l—以风景画的形式运行GrADS,此时其硬拷贝输出的区域大小为11×8.5英寸。 p—以肖像画的形式运行GrADS,此时其硬拷贝输出的区域大小为8.5 ×11英寸。 c—在GrADS启动后，首先执行其后提供的命令。 如果直接键入回车，GrADS将以风景画的形式启动。 *因此，在画图时要注意，不要把画图区域取得超过硬拷贝输出区域的大小。 *进入GrADS后一般可以看到两个窗口，上面一个为字符窗口，可以输入命令以及显示执行命令后的回应信息；下面一个为图形输出窗口。两个窗口中红色的为当前窗口。 *退出GrADS时，在字符窗口中键入：
一直听到有人抱怨，GrADS的坐标轴怎么那么固定，要设置个起始间隔还不支持时间轴，要在经纬度后面加个小圆圈的度数怎么就那么难，为啥不能四周都弄上坐标轴？ 好了，你的这些抱怨将会随着你看到这个帖子而消失，下面给出这个脚本的用法介绍： 注意：在你display变量前需要先设置如下两个命令 'set xlpos -20' 'set ylpos -20' 这样能够屏蔽GrADS自带的坐标轴，否则会出现坐标轴重叠 [code=gs] 请先display变量，再运行该脚本 该脚本主要实现了自定义的坐标轴显示，您可以免费使用该脚本 该脚本的用法规则如下： #axis -param1 val1 -param2 val2... 其中，-param表示需要自定义的参数类型，val表示该参数的具体设定值 参数和参数值必须成对出现，例如 #axis -type b -position o -sinterval 2 上面的语句表示绘制类型是底部的x轴，刻度位置是朝外，每两个大的刻度之间显示两个小刻度线 #注意：-type是必选参数 所有参数如下所示： -type ：表示绘制类型，参数值有：b/l/t/r 分别表示底部，左侧，上部，右侧 -label：表示是否绘制数值标签，参数值有：on/off 分别表示绘制和不绘制 -position：表示刻度的位置，参数值有：o/i/m 分别表示外侧，内侧，中间 -start：表示刻度从该处开始绘制，参数值为该轴范围内的数字 当该轴为时间轴时，表示开始绘制的时次（1,2,3...） -end：表示刻度的结束位置，参数值范围同start，时间轴时表示结束的时刻(1,2,3...) -interval：表示刻度的间隔，隔多少值绘制一个刻度和标签 例如，当x轴表示经度110-150E，如果interval为10，默认情况下会标出 110,120,130,140,150 不给定该参数时，默认情况下，会绘制5个刻度，以此确定默认间隔 -sinterval：表示没两个大刻度之间要绘制几个小的刻度，不给定该参数则不绘制 -size：表示刻度的长度，单位是英寸，默认为0.1 -color：表示坐标轴和刻度标签的颜色，默认为1 -lfont：表示坐标轴标签的字体，默认为当前环境下设置的字体 -lsize：表示坐标轴标签的字体大小，默认为0.12 -lthick：表示坐标轴标签字体的粗细，默认为0.3 -langle：表示坐标轴标签旋转的角度，默认为不旋转 -suffix：表示需要在每个坐标轴标签后面添加的文字或者符号后缀，比如用来添加度数等-asuffix：表示当坐标轴为经纬度轴时，是否自动添加E/W/N/S/EQ这种标记，默认为添加，参数值为 on/off ，分别表示自动添加和关闭自动添加 -hoffset：表示坐标轴标签在水平方向上的偏移量，单位为英寸，正负均可 -voffset：表示坐标轴标签在垂直方向上的偏移量，单位为英寸，正负均可 -tformat：当所绘制坐标轴为时间轴时产生作用，表示要显示的时间格式
以下技巧总结都是笔者从学习实践过程中总结出来的，基本的问题。不求全面，希望对读者学习有用，如果有问题，敬请留言指正，以促进交流学习！ 1、软件综述：grads软件是一款绘图软件除了绘制图形，还可以提取数据，主 要应用是在大气科学中，当然只要是数据处理成grads能够读取的数据文件就可以进行相关绘图。软件版本问题，软件本身不是很大，我接触到1.8、1.9、 2.0版本的，1.8版本的安装很多情况还要修改环境变量、1.9版本的不识别 ‘sdfopen’命令，最稳定的版本是2.0版本，所以笔者推荐学习者安装2.0版本，选择默认安装路径就可以。 2、文件类型简述：grads处理的是网格数据，可以处理的数据类型有：grd、 grib、nc（海洋常用的数据），cdf（雷达卫星数据），其中nc、cdf数据都是自带描述文件，不需要ctl，grib数据要通过命令生成ctl、index数据才可以调用，常用的是grd数据，需要ctl。 3、数据文件转换：grads软件识别的数据是二进制无格式数据，文件类型是 ‘binary’，写入和生成时是不需要格式的如read(20) sst(i,j,iz,it),20为文件号，通常是十进制数据与grd数据间转换，这里给一个grd转换成txt数据的fortran程序： parameter(nx=56,ny=41,nz=1,nt=360) dimension sst(nx,ny,nz,nt) real sst open(15,file='sst.grd',form='binary') ！固定的用form=‘binary’就是二进制数据open(16,file='sst.txt') ！新建txt文件 do it=1,nt do iz=1,nz read(15) ((sst(i,j,iz,it),i=1,nx),j=1,ny) ！read后只有文件号，数据是无格式的 enddo enddo do it=1,nt do iz=1,nz write(16,*) ((sst(i,j,iz,it),i=1,nx),j=1,ny) ！输出时是txt文件可直接看的数据，有格式输出，有* enddo enddo
dset ^test01.grb 给定二进制原始数据集的文件名为test01.grb。 title "Part I: test 01" 用字符串Part I: test 01 简略描述数据集内容 undef 1e+20 定义缺测值或缺值值，GrADS 在运算操作和图形操作时将忽略这些值点。dtype grib 数据集的数据类型为grib。 index ^test01.gmp 索引test01.gmp这个文件 xdef 72 linear 0.000000 5.000000 设置网格点值与经度的对应关系，x方向有72个网格点数 ydef 46 linear -90.000000 4.000000 定义y 方向格点与纬度的映象关系，表明共有46 个y 方向网格点 zdef 7 levels 设置气压面与垂直网格点的映射关系，共7层等压面，分别为1000hpa、850hpa、700hpa、500hpa、300hpa、200hpa、100hpa 1000 850 700 500 300 200 100 tdef 5 linear 0Z2jan2011 1dy设置网格值与时间的映射关系，表示共有5个时次，起始时刻2011 年1月2 号0Z 时，增量为1天。 vars 3 表示变量描述开始，有三个变量 q 7 51,100 Specific humidity [kg/kg] 第一个变量为q，有7层 ts 0 11,105, 2 Surface (2m) air temperature [K] 第二个变量为ts，有1层p 0 59, 1, 0, 0 Total precipitation rate [kg/(m^2*s)] 第三个变量为p，有1层endvars 表示数据描述文件结束
相关附件：(共323223 字节) funny给你一个小程序，是用perl写的，然后funny转成了exe文件，你可以用它生成ctl，但生成的ctl文件还需要自己去掉这个程序强制添加上去的一些信息，然后，你用gribmap.exe生成index文件，就可以显示了。 这个zip包里有原始的perl程序，转好的exe文件和gribmap.exe三个文件。 C:\drawing\ncep-monthly\ex>grib2ctl -i prs.grib.mean.y1980>y1980.ctl Using NCEP reanalysis table, see -ncep_opn, -ncep_rean options Using NCEP reanalysis table, see -ncep_opn, -ncep_rean options C:\drawing\ncep-monthly\ex>gribmap -i y1980.ctl Open Error: Unknown keyword in description file --> The invalid description file record is: --> this exe file was created with the evaluation version of perl2exe. The data file was not opened. File name is: y1980.ctl "Using NCEP reanalysis table, see -ncep_opn, -ncep_rean options" 就是说你应该用-ncep_opn 或-ncep_rean 的选项，具体看帮助 下面就不用说了，ctl都没有形成，自然不行的了 Hi,funny Thanks！ 错误与“-ncep_opn, -ncep_rean options”无关，是ctl中endvars后“--> this exe file was created with the evaluation version of perl2exe.”的这句话作怪，删掉后可正常得到idx。 ＝＝＝＝＝＝＝＝＝ D:\PCGrADS\win32>grib2ctl E:\。。。\air2m.mon.mean.nc >air2m.ctl 显示的信息如下： Big problem:
GrADS绘图软件安装于入门
目录 第一章GrADS绘图软件概述 1.GrADS绘图软件简介 2.GrADS绘图软件的安装（windows环境） 3.1在windows环境下安装GrADS软件包第二章GrADS绘图模板 1.GrADS示例演示 1.1 启动GrADS 1.2 退出GrADS 1.3 示例演示GrADS命令的使用
第二章GrADS绘图软件概述 1GrADS绘图软件简介 The Grid Analysis and Display System(GrADS) 是一套应用广泛、使用方便的科学数据绘图软件包。其主要特点： ●GrADS属于自由软件，可以从Internet上免费获得。 ●可运行于各种Windows 和Unix工作平台。 ●GrADS可用于4D数据的分析。既经度、纬度、层（气压层、高度层等） 和时间/xyzt 4维。数据可以是格点化的数据或离散点数据。GrADS 特别适用于气象类数据的分析。但也完全可以用于更广泛类型的数据分 析。 ●GrADS有多种显示方式：等值线、流线、矢量图、风矢量图、站点填 图、折线图、直方图等多种两维图形。 ●可处理多种数据格式的数据。GRIB、NetCDF、HDF-SDS等通用数据格 式和系统自定义的一种二进制数据格式。 ●采用命令行输入的方式交互式地显示图形。并有多种命令对数据进行再 加工。如求平均；计算涡度、散度、垂直积分、计算差分等。 ●图形可以按多种格式存储：ps、png、jpg、tiff、gif、windows metafile 等。自身以metafile格式存储。 ●提供多种辅助工具软件。如看图、打印、图形格式转换（gv）等工具。 2Internet上的GrADS资源 2.1 GrADS在Internet上的主页 GrADS主页地址： 从GrADS主页上可以找到预编译好了的适合于windows环境下的GrADS软件包， 2.2 windows环境下GrADS资源 ●grads-2.0.a7.oga.3-win32_superpack.exe——GrADS软件包安装执 行程序。目前在windows环境下较新的版本为GrADS第2.0版。 ●下载GrADS演示数据： 从下载： model.le.dat 和model.le.ctl文件； 或者GRIB码格式的数据： model.grb、model.gmp和model.ctl 3GrADS绘图软件的安装（windows环境） 在windows下运行GrADS的条件： ●安装GrADS软件包 3.1在windows环境下安装GrADS软件包 运行上述可执行文件（双击）(grads-2.0.a7.oga.3-win32_superpack.exe)进
GRADS (Grid Analysis and Display System) 讲义（6） 一、GRADS所能处理的数据格式 （1）无格式格式文件(直接、顺序存取)，Fortran可以读写。*.grd （2）N etCDF格式资料文件：精确性好，便于传输；*.nc （3）GRIB文件格式：压缩率高。 二、下面介绍第一种格式 GrADS中数据文件和数据描述文件是分开的。数据文件的存放一般为二进制直接访问形式(binary direct access)，其格式说明由数据描述文件(*.ctl)描述，该文件为纯文本格式，可用一般的编辑器产生（如EDIT,PE2等）。在GrADS环境中至少得首先打开(open)一个数据描述文件，以便后续的操作有数据对象。 1、首先搞清楚GrADS中数据的存放形式(五维的数据文件) （x，y）,z,VAR,t
GrADS格点为直接访问形式，一个网格点上(即一个确定的经纬度、高度和时刻)可以有任意多个物理变量，GrADS 视这些数据为一个大数组，其排放顺序为先经度、纬度、高度，然后是物理变量，最后是时次变化。一个x、y数据场构成一个记录，其顺序是x从西变到东，y从南变到北，从下到上，即实际大数组以二维数据片存放。 2、如何生成一个GrADS软件使用数据格式（sy.for） 现有ASCII码数据资料文件u.dat、v.dat和SST.dat，其空间范围60~150o E，0~40o N；层次：u、v为850、 200hPa；时段：1982.1~1985.12；分辨率：2.5*2.5。要求编写出将这三个文件转换成二进制无格式 直接存取（Grads格式）文件的Fortran程序，并给出相应的数据描述文件(CTL文件)。 C 定义一个X,Y,Z方向的格点以及总时次 nt parameter(nx=37,ny=17,nz=2,nt=48)
GrADS站点资料的使用 台站型离散资料画图的一般原则 1、把台站资料r.dat写成二进制文件r.grd 一般而言，台站资料都是文本格式的，须用专门的程序写成带有站号、经度、纬度等的二进制文件。 例如：对某一时次的降水资料r.dat有如下形式： lon lat Precipitation 119.8 30.6 0 119.8 30.616 0 119.8 30.632 0 119.8 30.648 0 119.8 30.664 0 119.8 30.68 0 119.8 30.696 0 119.8 30.712 0 119.8 30.728 0 119.8 30.744 0 119.8 30.76 0 119.8 30.776 0 。。。。。。 该数据只有一个时次 则相应的程序stn.f如下： parameter(n=10201) integer r1(n) real lat(n),lon(n),r(n) character*8 zh(n) open(1,file='d:\common\1.txt',form='formatted',status='old') do i=1,n read(1,*)lon(i),lat(i),r1(i) r(i)=r1(i) enddo do i=1,n zh(i)=char(i) enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!! open(9,file='d:\common\r.grd',form='binary',status='unknown') open(99,file='d:\common\rr.txt') do 100 i=1,n tim=0.0 nlev=1 nflag=1
G r ADS (Grid Analysis and Display System) 讲义（4） 内容提要 一、GrADS所能处理的数据格式 二、数据资料的准备（binary） ASCII码数据文件向二进制数据文件的转换 三、数据描述文件的构成 四、其它数据文件简介（netCDF/GRIB） 五、GrADS控制文件简介
一、GrADS所能处理的数据格式 --Binary：自制，直接、顺序存取，Fortran 可以读写。*.grd --netCDF（自描述）：精确性好，便于传输；*.nc --GRIB：压缩率高。*.grb --HDF-SDS（卫星资料） --ASCII（台风路径） --站点 --BUFR（V1.9版本）：二进制通用数据表示格式（BUFR），用于非格点气象数据的保存，便于网络传输，是世界气象组织（WMO）规定的标准格式，目前使用的常规气象资料数据以及雷达、卫星数据资料转换为BUFR格式 二、Binary格式介绍 1、简介 1）GrADS最基本的数据格式； 2）常用后缀：*.grd,.bin,.dat；
3）可用fortran读写； 4）一般为二进制无格式文件(form=‘unformatted’)； 5）访问形式 直接（access=‘direct’） 顺序（access=sequential）--ctl文件中说明6）格式说明由数据描述文件(*.ctl)描述； 数据资料———数据描述文件———GrADS 翻译器 注意：数据文件和数据描述文件是分开的。后者为纯文本格式，可用一般的编辑器产生（如记事本等）； 2、GrADS中数据的存放形式 1）5-D的数据集 （x，y）,z,VAR,t
GrADS绘图软件实用手册 2002年1月
目录 第一章GrADS绘图软件概述 1.GrADS绘图软件简介 2.Internet上的GrADS资源 2.1GrADS在Internet上的主页 2.2windows环境下GrADS资源 3.GrADS绘图软件的安装（windows环境） 3.1在windows环境下安装GrADS软件包 X server 的安装 第二章GrADS绘图模板 1.GrADS示例演示 1.1 启动GrADS 1.2 退出GrADS 1.3 示例演示GrADS命令的使用 2.GrADS绘图模板 3.GrADS模板的高级应用 GrADS描述语言 GrADS高级模板的应用 第三章GrADS数据格式 1.格点数据描述文件 1.1 数据描述文件各项解释 1.2 生成model.le.dat和model.le.ctl文件的程序代码片段 2.站点数据的格式 附录 1.如何精确控制图形输出的尺寸—Landscape纸型 2.台站资料的显示 3.Linux环境下的安装
第二章GrADS绘图软件概述 1GrADS绘图软件简介 The Grid Analysis and Display System(GrADS) 是一套应用广泛、使用方便的科学数据绘图软件包。其主要特点： ●GrADS属于自由软件，可以从Internet上免费获得。 ●可运行于各种Windows 和Unix工作平台。 ●GrADS可用于4D数据的分析。既经度、纬度、层（气压层、高度层等） 和时间/xyzt 4维。数据可以是格点化的数据或离散点数据。GrADS 特别适用于气象类数据的分析。但也完全可以用于更广泛类型的数据分 析。 ●GrADS有多种显示方式：等值线、流线、矢量图、风矢量图、站点填 图、折线图、直方图等多种两维图形。 ●可处理多种数据格式的数据。GRIB、NetCDF、HDF-SDS等通用数据格 式和系统自定义的一种二进制数据格式。 ●采用命令行输入的方式交互式地显示图形。并有多种命令对数据进行再 加工。如求平均；计算涡度、散度、垂直积分、计算差分等。 ●图形可以按多种格式存储：ps、png、jpg、tiff、gif、windows metafile 等。自身以metafile格式存储。 ●提供多种辅助工具软件。如看图、打印、图形格式转换（gv）等工具。2Internet上的GrADS资源 2.1 GrADS在Internet上的主页 GrADS主页地址：https://www.doczj.com/doc/9c12703299.html,/grads 从GrADS主页上可以找到预编译好了的适合于windows环境下的 GrADS软件包，和适合于各种UNIX环境下的GrADS软件包。关于 GrADS在UNIX环境下的资源请参考附录。 2.2 windows环境下GrADS资源 ●GrADS1.8sl8.win32.exe——GrADS软件包安装执行程序。目前在 windows环境下最新的版本为GrADS第1.8版。 ●下载GrADS演示数据： 从ftp://https://www.doczj.com/doc/9c12703299.html,/grads/sprite/tutorial下载： model.le.dat 和model.le.ctl文件； 或者GRIB码格式的数据： model.grb、model.gmp和model.ctl 3GrADS绘图软件的安装（windows环境） 在windows下运行GrADS的条件： ●安装GrADS软件包 ●安装X SERVER软件包 3.1在windows环境下安装GrADS软件包 运行上述可执行文件(GrADS1.8sl8.win32.exe)进入第一个画面：