台风路径频数绘制 MATLAB+GrADS

dieam

之前画了多年的台风路径  ----->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. 
    
24. for i=1:282  %台风数************************************************
    
25.     ii=1;
    
26.     for j=1:144 %格点数 144*25=3600
    
27.         jj=1;
    
28.        for k=1:36 %格点数 36*25=900
    
29.            D(j,k,i)=sum(sum(C(ii:j*25,jj:k*25,i)));
    
30.            jj=jj+25;
    
31.        end
    
32.        ii=ii+25;
    
33.     end
    
34. end
    
35. sum(sum(sum(D(:,:,:))))
    
36. for i=1:282  %台风数**********************************************
    
37.     for j=1:144
    
38.         for k=1:36
    
39.             if D(j,k,i)~=0
    
40.                 DD(j,k,i)=1;
    
41.             else
    
42.                 DD(j,k,i)=0;
    
43.             end
    
44.         end
    
45.     end
    
46. end
    
47. for i=1:144
    
48.     for j=1:36
    
49.         DDD(i,j)=0;
    
50.         for k=1:282 %台风数***********************************************
    
51.             DDD(i,j)=DD(i,j,k)+DDD(i,j);
    
52.         end
    
53.     end
    
54. end
    
55. x=40:1:100;
    
56. y=1:1:28;
    
57. contour(x,y,DDD(40:100,1:28)')
    
58. 
    
59. %A=xlsread('d:/bylwdata/png/yeonumber.xlsx','Sheet1');
    
60. %  AA=flipud(DDD);
    
61. %  AA=AA';
    
62. fid=fopen('E:/grads/data/tcpath.dat','wb');
    
63. count=fwrite(fid,DDD,'float32');
    
64. 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出的图

grads平滑前

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

我也把三个文件打包上传
path program.zip (1.42 KB, 下载次数: 22)

2020-4-2 22:07 上传

点击文件名下载附件

伽蓝鸟

我之前写过一个脚本，让MATLAB读入热带气旋最佳路径集的数据，你在二爷的教程里能找到。另外，你知道有个函数叫histogram2咩。。

dieam

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

我用matlab比较少，刚看到这个函数有些懵，所以才和同学讨论写出来这个很麻烦的方法。
我现在在看您写的脚本，也在看这个函数的使用方法，争取能写出来一个比较好的程序

climate777

同学你好，我感觉我们毕业论文方向好像，我可以加你qq吗？有些问题想请教一下

997855rr

哇 感谢楼主 感觉很有用

小谢

hello，请问可以加一下QQ或者微信么？
想请教一下关于台风路径频数的问题

miceyou

非常感谢分享~！太有用了！

Ruby_Guuu

感谢楼主分享~