GrADS台风路径gs脚本分享，个人修正版（资源+样图）

ljh110011

GrADS台风路径搜寻作图资源分享

by ljh110011

2015.05.23

这个grads脚本文件是参考网上一些台风搜寻路径一些资源，自己作很多更新和修正的。增加了很多个性化的标注，图样等等。

非常感谢网上各前辈的资源和给予的参考。

下面贴出样图和脚本代码。脚本代码lj1.gs和lj2.gs，还有测试数据TestData.tar.gz附件中也可以下载得到。

希望各位觉得本帖方法还可以的，多多回帖支持和指正错误，谢谢。

个人写的WRF数值模拟模式系列软件安装教程链接：

http://bbs.06climate.com/forum.php?mod=viewthread&tid=33557

提供一个好地方：

中国气象局热带气旋资料中心：

可以下载最佳台风路径数据集和卫星分析资料等很多资料。

http://tcdata.typhoon.gov.cn/index.html

万能胶卷

有理有据，非常清晰

ljh110011

下面是lj1.gs脚本的命令：

*从海平面气压场（或1000hPa高度场）中搜寻台风中心位置，输出到文件

'reinit'

'!rm moni-w.txt'  ;*使用追加写入，所以要先删除旧数据文件 !是GRADS调用操作系统命令

'open /home/DATA/a/lj1.ctl'   ;*打开经过ARWpost转换后WRF预报结果 只需要slp变量

i=1

while (i<8) ;*要搜索的dat文件中n个时刻的中心点 填写n+1

'set t ' i

t=subwrd(result,4) ;*取得日期 result截取数据中第4个单词

tc=substr(t,6,7) ;*substr函数 变量 起始字符 字符长度

'd slp'  ;*处理海平面气压变量

'definelocx=minloc(min(slp,lat=10,lat=20),lon=115,lon=140)' ;*搜索最低气压经度 可调节

'definelocy=minloc(min(slp,lon=115,lon=140),lat=10,lat=20)' ;*搜索最低气压纬度 可调节

'd locx'

locxx=subwrd(result,4)

'd locy'

locyy=subwrd(result,4)

'q gr2w 'locxx' 'locyy ;*GRADS坐标转换成经纬度坐标

x0=subwrd(result,3)

y0=subwrd(result,6)

*控制输出数据有几个小数位 实在麻烦 其实可以用excel作数据处理

*保留一位小数设置

x=math_format('%3.1f',x0)
y=math_format('%2.1f',y0)

rc=write('moni-w.txt',x ' 'y '  ' tc,append)  ;*将经纬度位置和时间存入文件

i=i+1   ;*这里灵活运用 根据dat文件是多少小时的数据设置

endwhile

close('moni-w.txt') ;*关闭输出文件

'!cp moni-w.txt moni.txt' ;*在GS里同一个文件不能既读又写，所以做1个拷贝

下面是lj2.gs脚本的命令：

'reinit' ;*从上述数据文件读出位置绘制台风路径

'open /home/DATA/a/lj1.ctl'

'set grads off'

'set mproj latlon'

'set poli on'     ;*画出国界省界线

'set mpdset hires' ;*设置高分辨率地图

*定义设置排版要求

'set lon 105 145'  ;*太宽纵坐标显示缺部分

'set lat 0 30'

'set cmax 0'  ;*把显示最大值设小，是为了只出底图

'set xlint 5'  ;*定义横坐标的标记间隔

'set ylint 5'  ;*定义纵坐标的标记间隔

'set xlopts 1 4 0.15' ;*横坐标 颜色 线宽 大小

'set ylopts 1 4 0.15' ;*纵坐标 颜色 线宽 大小

'set clopts -1 -1 0.08' ;*等值线 颜色 线宽 标记大小

'd slp' ;*作出底图用

*开始画出模拟数据路径

'set string 2 bl 1 60'   ;*设置字符串颜色  位置  粗细  旋转角度

'set strsiz 0.06 0.06'    ;*时间标注的线宽  大小

'set line 2 1 0.1'    ;*设置线条类型红色 实线 粗细

*使用以上定义，基本符合排版图约5*7CM大小时图中标注约为6号字大小

i=1

while(i<=15)  ;*文件内容n行之内

read_file=read('moni.txt')

read_code=sublin(read_file,1)  ;*读取文件打开代码

if(read_code>0);

break;

endif ;*读取正常为0，其他要么出错，要么结束

xbefore=x0

ybefore=y0

read_line=sublin(read_file,2) ;*读取文件行内容

lon0=subwrd(read_line,1)  ;*取得第1部分 经度

lat0=subwrd(read_line,2)  ;*取得第2部分 纬度

date0=subwrd(read_line,3)  ;*取得第3部分 时间

'q w2xy ' lon0 ' ' lat0  ;*将经纬度坐标转换为GRADS坐标

x0=subwrd(result,3)  ;*取得转换后的X坐标

y0=subwrd(result,6)  ;*取得转换后的Y坐标

'draw mark 3 'x0' 'y0' 0.1'   ;*作出各个坐标点

if(i>1) ;*第一个点作为起始点

'draw line 'xbefore' 'ybefore' 'x0''y0''  ;

endif ;

'draw string 'x0+0.1' 'y0+0.1' 'date0'('lat0','lon0')'  ;*标注点的坐标和信息

*if(math_mod(i,2));'draw string 'x0-0.1''y0+0.2' 'date0'';endif;

i=i+1

endwhile

*画出图样和注释

'draw mark 3 1 0.3 0.1'  ;*第一个点

'draw mark 3 1.5 0.3 0.1'  ;*第二个点

'draw line 1 0.3 1.5 0.3'  ;*两点连线

'set font 2'    ;*选择字形库  0到5

'set string 2 c 2 0'   ;*设置字符串颜色  位置  粗细  旋转角度

'set strsiz 0.11  0.15'  ;*设置字符串的水平大小和高度大小

'draw string 2.25 0.3 Analog track'   ;*设置脚下标

*开始画出实况数据路径

'set string 4 br 1 60'   ;*设置字符串颜色  位置  粗细  旋转角度

'set strsiz 0.06 0.06'    ;*时间标注的线宽  大小

'set line 4 3 0.1'    ;*设置线条类型蓝色 短虚线 粗细

*使用以上定义，基本符合排版图约5*7CM大小时图中标注约为6号字大小

i=1

while(i<=15)  ;*文件内容n行之内

read_file=read('shikuang.txt')

read_code=sublin(read_file,1)  ;*读取文件打开代码

if(read_code>0);

break;

endif ;*读取正常为0，其他要么出错，要么结束

xbefore=x0

ybefore=y0

read_line=sublin(read_file,2) ;*读取文件行内容

lon0=subwrd(read_line,1)  ;*取得第1部分 经度

lat0=subwrd(read_line,2)  ;*取得第2部分 纬度

date0=subwrd(read_line,3)  ;*取得第3部分 时间

'q w2xy ' lon0 ' ' lat0  ;*将经纬度坐标转换为GRADS坐标

x0=subwrd(result,3)  ;*取得转换后的X坐标

y0=subwrd(result,6)  ;*取得转换后的Y坐标

'draw mark 3 'x0' 'y0' 0.1'   ;*作出各个坐标点

if(i>1) ;*第一个点作为起始点

'draw line 'xbefore' 'ybefore' 'x0''y0''  ;

endif ;

'draw string 'x0-0.1' 'y0-0.1' 'date0'('lat0','lon0')'  ;*标注点的坐标和信息

*if(math_mod(i,2));'draw string 'x0-0.1''y0+0.2' 'date0'';endif;

i=i+1

endwhile

*画出图样和注释

'draw mark 3 5 0.3 0.1'  ;*第一个点

'draw mark 3 5.5 0.3 0.1'  ;*第二个点

'draw line 5 0.3 5.5 0.3'  ;*两点连线

'set font 2'    ;*选择字形库  0到5

'set string 4 c 2 0'   ;*设置字符串颜色  位置  粗细  旋转角度

'set strsiz 0.11  0.15'  ;*设置字符串的水平大小和高度大小

'draw string 6.25 0.3 Sence track'   ;*设置脚下标

*设置标题,标注等

'set string 1 c 2 0'   ;*设置字符串颜色  位置  粗细  旋转角度

'set font 2'    ;*选择字形库  0到5

'set strsiz 0.11  0.15'  ;*设置字符串的水平大小和高度大小

'draw string 9.8 0.3 2015-05-23ljh110011'   ;*设置脚下标

'set strsiz 0.3  0.6'  ;*设置字符串的水平大小和高度大小

'draw title 2013-08-09~2013-08-16\TyphoonUtor Track'    ;*在图片顶部写字符串作为标题

'printim /home/DATA/a/l3.png white'  ;*图形输出

至此，gs脚本程序贴出完毕，希望各位多多支持和指正，谢谢。

river

感谢分享，很好的脚本，标注也很详细

兰溪之水

控制格式化输出，有相应的函数，如
LonM = math_format('%5.1f',LonM)

Jepeng

好厉害，谢谢楼主分享

任我行

thank you!!!!!!!!!!!!!!!!!!!!!!!!!!

ljh110011

兰溪之水 发表于 2015-5-23 19:52
控制格式化输出，有相应的函数，如
LonM = math_format('%5.1f',LonM)

谢谢兰溪大大指导，还以为grads不能自动取小数位，现在学会了。

ljh110011

kakame08 发表于 2015-5-24 15:44
这个很棒啊
是不是可以画任意坐标路径？
一直很需要

可以的，把经纬度列在txt文件下，然后引到lj2.gs的路径下就能实现了。

zhaolanxiai

牛人啊！其实如果有路径数据，用surfer也可以做的