GRADS绘制台风(热带气旋)路径-自动识别中心-位置和最低气压保存

蕎蕎安

{:5_213:}{:5_213:}

gingili

执行之后图不对，图片不对标注也不对
'reinit' ;*第1部分--从海平面气压场（或1000hPa高度场）中找出台风中心位置，输出到文件
'!rm weizhi-w.txt' ;*使用追加写入，所以要先删除旧数据文件
'open E:\20160503\datafiles1\fnl_20160501_00_00.ctl';*WRF预报结果
i=1
'set lev 1000'
while (i<16) ;*本来应该用<=,但是最后一个位置不对，不用
'set t 'i
t=subwrd(result,4) ;*取得日期
tc=substr(t,8,5) ;*取得部分日期串
'd HGTprs';*海平面气压
'define locx=minloc(min(HGTprs,lat=30,lat=55),lon=110,lon=140)' ;*搜索最低气压经度
'define locy=minloc(min(HGTprs,lon=110,lon=140),lat=30,lat=55)' ;*搜索最低气压纬度
'd locx'
locxx=subwrd(result,4)
'd locy'
locyy=subwrd(result,4)
'q gr2w 'locxx' 'locyy ;*GRADS坐标转换成经纬度坐标
x=subwrd(result,3)
y=subwrd(result,6)
'set x 'locxx
'set y 'locyy
'd HGTprs'
minslp=subwrd(result,4)
rc=write('weizhi-w.txt',x ' 'y ' 'minslp ' 'tc,append)  ;*将经纬度位置、最低气压和时间存入文件
*say rc
i=i+4 ;*只输出08和20时，数据间隔3小时
endwhile
close('weizhi-w.txt') ;*关闭输出文件
'!cp weizhi-w.txt weizhi-r.txt' ;*在GS里同一个文件不能既读又写，所以做1个拷贝

'reinit' ;*第2部分--从上述数据文件(你也可自己编辑)读出位置绘制台风路径
'open E:\20160503\datafiles1\fnl_20160501_00_00.ctl'
'set mproj latlon'
'set mpdset hires' ;*设置高分辨率地图
'set lon 110 140' ;*太宽纵坐标显示缺部分
'set lat 30 55'
'set lev 1000'
'set cmax 0' ;*把显示最大值设小，是为了只出底图

'set xlint 2' ;*定义横坐标的标记间隔
'set ylint 2' ;*定义纵坐标的标记间隔
'set xlopts 1 4 0.22' ;*横坐标 颜色 线宽 大小
'set ylopts 1 4 0.20' ;*纵坐标 颜色 线宽 大小
'set clopts -1 -1 0.20' ;*等值线 颜色 线宽 标记大小
*使用以上定义，基本符合排版要求--图大约5*7CM坐标轴标注约为6号字大小

'd HGTprs' ;*只为出底图用
read_file=read('weizhi-r.txt') ;*打开文件进行读取
read_line=sublin(read_file,2) ;*返回结果分为两行（但仍作为一个串）
lon0=subwrd(read_line,1) ;*读取第1行中经度
lat0=subwrd(read_line,2) ;*读取第1行中纬度
date0=subwrd(read_line,3) ;*读取第1行中时间
'q w2xy ' lon0 ' ' lat0 ;*将经纬度转换成GRADS坐标
x0=subwrd(result,3) ;*取得转换后的X坐标
y0=subwrd(result,6) ;*取得转换后的Y坐标
'draw mark 3 'x0' 'y0' 0.12' ;*先标注起始点
'set string 2 '
'set strsiz 0.18 0.22' ;*时间标注的 线宽 大小
'draw string 'x0-0.1' 'y0+0.2' 'date0''
*使用以上定义，基本符合排版图约5*7CM大小时图中标注约为6号字大小

i=1
while(i<=200) ;*文件内容200行之内
xbefore=x0
ybefore=y0
read_file=read('weizhi-r.txt')
read_code=sublin(read_file,1) ;*读取文件打开代码
if(read_code>0);break;endif;*读取正常为0，其他要么出错，要么结束
read_line=sublin(read_file,2) ;*读取文件行内容
lon0=subwrd(read_line,1) ;*取得第1部分
lat0=subwrd(read_line,2) ;*取得第2部分
date0=subwrd(read_line,4) ;*取得第3部分
'q w2xy ' lon0 ' ' lat0 ;*将经纬度坐标转换为GRADS坐标
x0=subwrd(result,3) ;*取得X坐标
y0=subwrd(result,6) ;*取得Y坐标
*'draw title 1119 "nige" track '
'draw mark 3 'x0' 'y0' 0.12'
'draw line 'xbefore' 'ybefore' 'x0' 'y0''
*if(math_mod(i,2));'draw string 'x0-0.1' 'y0+0.2' 'date0'';endif;
'draw string 'x0-0.1' 'y0+0.2' 'date0'';
i=i+1
endwhile
'gxprint track.png white';*图形输出
;

                               
登录/注册后可看大图

幽居寒舍

牛逼Le，学习

gingili

楼主，执行程序之后中间显示unknow command：0，注释掉close('weizhi-w.txt') ;*关闭输出文件，就没有这个问题，想问一下影响最终结果不？

xxmoumou

看着代码一脸懵，要好好学学呢！

熙熙攘攘

谢谢提供，正好想看看

NerryGIS

{:eb513:}{:eb513:}{:eb513:}{:eb513:}{:eb513:}{:eb513:}{:eb513:}

jsl

学习一下，感谢

zhyychbyzh

还没看到，先谢谢了

小明同学

感谢楼主啦！！！