请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 20826|回复: 7

[分享资料] 台风路径频数绘制 MATLAB+GrADS

[复制链接]

新浪微博达人勋

发表于 2020-4-2 22:01:09 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 dieam 于 2020-4-2 22:07 编辑

之前画了多年的台风路径  ----->grads画多条台风路径   虽然能够用不同的颜色表示强度,但最后画出来还是太过密集,不容易分析。和老师交流之后,老师让我画台风路径频数的等值线图。即:经纬度划分成2.5°×2.5°的格点,路径每经过一次就+1,统计多年路径频数。
之前在家园搜索,有前辈用MeteoInfo进行绘制 http://bbs.06climate.com/forum.php?mod=viewthread&tid=68433&highlight=%C2%B7%BE%B6%C6%B5%CA%FD  ,但由于没有接触过MeteoInfo便放弃了。不过我也强烈推荐大家去学习。王老师的MeteoInfo功能很强大的!
-------------------------------------------------------------------------------------------
之后我也思考了好多天,终于有了结果,来把思路分享给大家。

我在绘制时,先用matlab处理了数据,然后输出成dat,编写ctl文件,最后通过grads进行绘制。

首先将将多年的台风路径数据读到excel里面,通过matlab进行计算。
在台风路径文件中,头记录中第6列为终结记录,0表示消散,所以通过判断此数据为0来识别不同的台风,if (B(i,3))==0。当然,由于在excel中分列时设置的不同,此项也可能为nan,此时要改为if isnan(B(i,3))。  要根据读入数据不同情况进行修改!!

路径数据中,经纬度的单位为0.1°,在编程时我也没有通过*0.1来换算,所以程序中经纬度设置3600和900(北半球)。相应的,2.5°×2.5°就变成了25×25。当然,这并不影响格点数,依旧为144和36。


整体的计算思路就是:如果读入的数据为A(i,j,k),K代表第三维,即要画的台风路径数量,一个台风路径数据存在一个二维数组里,经过那个格点,那个格点的数值就为1,没有经过的格点数据为0,然后利用sum函数进行区域频数统计,依旧得出A(i,j,k)的数据,但是i,j的范围变了,因为划分区域了,i,j变小了,然后这个时候进行对非0数值变为1,最后再把每一个台风数组的每个格点相加,得出要画的等值线的图的数据。
  1. clear
  2. clc
  3. A=xlsread('F:\MCM\data1.xlsx');
  4. B(:,1)=A(:,3); %读第3列  纬度
  5. B(:,2)=A(:,4); %读第4列  经度
  6. B(:,3)=A(:,6); %读第6列  在头记录里是0用来区别每个台风
  7. i=1;
  8. k=0;
  9. A1=max(B(:,1));
  10. A2=min(B(:,1));
  11. A3=max(B(:,2));
  12. A4=min(B(:,2));
  13. C(1:3600,1:900,1:282)=0; %经纬度和台风数**************************************
  14. for i=1:8716   %一共8716行数据 ***************************************
  15.     if (B(i,3))==0   
  16.         i=i+1;
  17.         k=k+1;
  18.     end
  19.     m=B(i,1);
  20.     n=B(i,2);
  21.     C(n,m,k)=1;
  22. end

  23. for i=1:282  %台风数************************************************
  24.     ii=1;
  25.     for j=1:144 %格点数 144*25=3600
  26.         jj=1;
  27.        for k=1:36 %格点数 36*25=900
  28.            D(j,k,i)=sum(sum(C(ii:j*25,jj:k*25,i)));
  29.            jj=jj+25;
  30.        end
  31.        ii=ii+25;
  32.     end
  33. end
  34. sum(sum(sum(D(:,:,:))))
  35. for i=1:282  %台风数**********************************************
  36.     for j=1:144
  37.         for k=1:36
  38.             if D(j,k,i)~=0
  39.                 DD(j,k,i)=1;
  40.             else
  41.                 DD(j,k,i)=0;
  42.             end
  43.         end
  44.     end
  45. end
  46. for i=1:144
  47.     for j=1:36
  48.         DDD(i,j)=0;
  49.         for k=1:282 %台风数***********************************************
  50.             DDD(i,j)=DD(i,j,k)+DDD(i,j);
  51.         end
  52.     end
  53. end
  54. x=40:1:100;
  55. y=1:1:28;
  56. contour(x,y,DDD(40:100,1:28)')

  57. %A=xlsread('d:/bylwdata/png/yeonumber.xlsx','Sheet1');
  58. %  AA=flipud(DDD);
  59. %  AA=AA';
  60. fid=fopen('E:/grads/data/tcpath.dat','wb');
  61. count=fwrite(fid,DDD,'float32');
  62. fclose(fid);
复制代码

最后输入成dat的时候要注意数据的格式,可以参考  matlab数据转grads格点数据函数  和  Matlab处理后的数据写为二进制用于GrADS画图——fwrite函数使用  两个帖子,不然画出来的图会有问题

我虽然通过这个方法计算出来路径频数,不过并不是最好的。txt读入excel,统计台风数量都很麻烦,而且在matlab中计算量也很大,如果统计的台风数量过多,会超出数组维数,无法计算。
可能histogram2计算会更方便,但我还没有写出来。 研究出来之后我也会及时分享给大家。

---------------------------------------------------------------------------------------------------------------------------------
输出dat之后就很方便了,编写一个ctl文件,就可以在grads中绘图了。只有一个变量,随便设置一个名字用于画图就好。lev和time我也都随意设置了一个。
  1. dset E:\data\tcpath1.dat
  2. title typhoon path
  3. undef -999
  4. xdef 144 linear 0 2.5
  5. ydef 36 linear 0 2.5
  6. zdef 1 levels 1000
  7. tdef 1 linear JAN2009 1yr
  8. vars 1
  9. path1 0 99
  10. endvars
复制代码
---------------------------------------------------------------------------------------------------------------------------------
在绘图时,我做了一个九点平滑,这样出来的图会更好看一点。
  1. 'reinit'
  2. 'open E:\data\tcpath1.ctl'
  3. 'set grid off'
  4. 'set grads off'
  5. 'set mpdset hires cnworld'
  6. 'set poli on'
  7. 'set map 1'
  8. 'set t 1'
  9. 'set z 1'
  10. 'set lon 90 170'
  11. 'set lat 0 55'
  12. 'set xlevs 90 100 110 120 130 140 150 160 170'
  13. 'set ylevs 0 5 10 15 20 25 30 35 40 45 50 55'
  14. 'set xlopts 1 0.1 0.16'
  15. 'set ylopts 1 0.1 0.16'
  16. 'set gxout contour'
  17. 'set cmin 1'
  18. 'set ccolor rainbow'
  19. 'set cthick 4'
  20. 'd smth9(path1)'
  21. 'draw title TC Pathway'
  22. 'gxprint E:\figure\path2.png white'
  23. ;
复制代码

------------------------------------------------------------------------------------------------------------------------------
最后,grads画出来的图可以和MATLAB出的图对比一下,看是否正确。

MATLAB出的图

MATLAB出的图

grads平滑前

grads平滑前
path2.png
可以看到,grads的图和MATLAB的图是一样的,所以我们输出数据及画图过程没有出问题。 不过九点平滑之后,应该更好看一点了。

我也把三个文件打包上传
path program.zip (1.42 KB, 下载次数: 21)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2020-4-3 08:11:07 | 显示全部楼层
我之前写过一个脚本,让MATLAB读入热带气旋最佳路径集的数据,你在二爷的教程里能找到。另外,你知道有个函数叫histogram2咩。。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-4-3 09:18:25 | 显示全部楼层
伽蓝鸟 发表于 2020-4-3 08:11
我之前写过一个脚本,让MATLAB读入热带气旋最佳路径集的数据,你在二爷的教程里能找到。另外,你知道有个函 ...

我用matlab比较少,刚看到这个函数有些懵,所以才和同学讨论写出来这个很麻烦的方法。
我现在在看您写的脚本,也在看这个函数的使用方法,争取能写出来一个比较好的程序
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-4-23 19:53:24 | 显示全部楼层
同学你好,我感觉我们毕业论文方向好像,我可以加你qq吗?有些问题想请教一下
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-4-13 19:17:52 | 显示全部楼层
哇 感谢楼主 感觉很有用
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-4-17 17:43:03 | 显示全部楼层
hello,请问可以加一下QQ或者微信么?
想请教一下关于台风路径频数的问题
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-4-18 09:07:22 | 显示全部楼层
非常感谢分享~!太有用了!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-4-14 16:39:47 | 显示全部楼层
感谢楼主分享~
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表