畫出風標時間序列圖 更新:上傳數據

river

传点儿数据，可以让大家帮你先试一下

ml11178

mofangbao 发表于 2012-10-15 10:44

                               
登录/注册后可看大图

额。。楼主是想用Fortran画还是你下面说的其他语言呢？Fortran画的话是比较痛苦的，不建议。
ncl我也不熟悉 ...

恩，沒錯，目前是希望先用FORTRAN先學著做，"言深深"大大提到的SUFER我也是很有興趣
我先上傳我的數據

麻煩各位大大指導小弟了

2012081201calwswd.dat (12.96 KB, 下载次数: 11)

2012-10-15 11:19 上传

点击文件名下载附件

mofangbao

传说中的谁 发表于 2012-10-15 10:47

                               
登录/注册后可看大图

貌似是上次你给我弄的那个玩意儿吧？
PS：欢迎宝岛的同行加入到论坛

恩 应该差不多 他这个是气球的  你那个时雷达的

mofangbao

ml11178 发表于 2012-10-15 11:20

                               
登录/注册后可看大图

恩，沒錯，目前是希望先用FORTRAN先學著做，"言深深"大大提到的SUFER我也是很有興趣
我先上傳我的數據
...

数据的每一列的含义介绍下呀

ml11178

mofangbao 发表于 2012-10-15 11:57

                               
登录/注册后可看大图

数据的每一列的含义介绍下呀

t是時間，每10秒一筆資料
VE、VN是兩個速度的分量
THETA是風向的角度
v為風速

mofangbao

ml11178 发表于 2012-10-15 12:04

                               
登录/注册后可看大图

t是時間，每10秒一筆資料
VE、VN是兩個速度的分量
THETA是風向的角度

如果需要显示随高度的分布的话，高度值应该要在数据里面的，你的数据没有高度值啊，比如t=0对应的高度是多少，每一个t都应该对应一个高度值，可以使高度值或者是气压值

mofangbao

你输出数据的时候 t=  vn=这种字符串就不用输出了，直接输出一列列的数据就行了，字符串和数据混合输出不利于Fortran再次读取

ml11178

mofangbao 发表于 2012-10-15 12:08

                               
登录/注册后可看大图

如果需要显示随高度的分布的话，高度值应该要在数据里面的，你的数据没有高度值啊，比如t=0对应的高度是多 ...

恩，因為只是簡單的想看風向的日夜變化，隨高度的變化不算是主要
所以沒有特別準確的算出高度或氣壓值

大概是下圖這樣，我猜下圖的風標是直接以每筆資料堆疊上去的，如果每筆資料的秒數一樣的話

高度大概也不會差多少，因為這個實驗的假設是氣球的垂直上升速度是固定的

縱軸為大概的高度，橫軸為小時

ml11178

mofangbao 发表于 2012-10-15 12:18

                               
登录/注册后可看大图

你输出数据的时候 t=  vn=这种字符串就不用输出了，直接输出一列列的数据就行了，字符串和数据混合输出不利 ...

恩，以後會多注意 :)

mofangbao

由于只有你的一个时次的数据，因此我把你的一个文件进行了复制，总共6个时次（当然数据都是一样的），来进行测试，思路如下：
先已经把你数据中的字符串删掉了，你重新生成的时候不要写入字符串了，处理后的数据如下：
   0.  -367.50  -889.42   67.55   96.24
  10.   357.95   776.84  245.26   85.53
  20.    21.46    30.93  235.24    3.76
  30.    29.86   -13.05  156.39    3.26
  40.    10.89   -26.41  112.41    2.86
  50.    14.47   -33.27  113.50    3.63
  60.    27.57   -27.52  135.05    3.90
  70.    29.75   -23.17  142.09    3.77
  80.    28.86   -47.56  121.25    5.56
  90.    27.77   -52.33  117.96    5.92
100.    24.57   -36.61  123.87    4.41
110.    25.74   -43.51  120.61    5.06
...
然后Fortran对数据文件进行提取和转换，生成GrADS格式的数据（可能你需要先学习一些GrADS的基础知识），只要使用到其中的最后两列，Fortran程序如下：

1. program wnd
   
2. !zcount表示一共高度上有多少层，你的数据后面有缺测，因此规定一个数值，保证每个文件读取到的层数为这个值，可自行修改
   
3. !undef为数据中没有这么多层的时候用来填补的缺测值
   
4. !startname和endname分别为数据开始和结束的时间，用来动态生成文件名，可自行修改
   
5. !程序执行后会分别生成theta和v的两个数据文件，每个文件中都有N个时次，注意该程序每次最多只能运行01时到23时，如果要更多时次，请自行修改程序中文件名的动态生成方法。
   
6. integer,parameter::zcount=100,undef=32766,startname=2012081201,endname=2012081206
   
7. real var(5)
   
8. character cname*21
   
9. integer::status=0,iname
   
10. 
    
11. open(2,file='20120812.theta.grd',status='replace',form='binary')
    
12. open(3,file='20120812.v.grd',status='replace',form='binary')
    
13. 
    
14. do iname=startname,endname
    
15. call getCname(iname,cname)
    
16. open(1,file='file\'//cname,status='old')
    
17. do i=1,zcount
    
18.   read(1,*,iostat=status)var
    
19.   if(status/=0)then
    
20.     write(2)undef
    
21.     write(3)undef
    
22.   else
    
23.     write(2)var(4)
    
24.     write(3)var(5)
    
25.   endif
    
26. enddo
    
27. close(1)
    
28. enddo
    
29. close(2)
    
30. close(3)
    
31. end
    
32. 
    
33. subroutine getCname(iname,cname)
    
34. integer::iname
    
35. character*21 cname
    
36. cname=''
    
37. write(cname,'(i10)')iname
    
38. cname=trim(cname)//'calwswd.dat'
    
39. return
    
40. endsubroutine
    
41. 
    

program wnd<br /> !zcount表示一共高度上有多少层，你的数据后面有缺测，因此规定一个数值，保证每个文件读取到的层数为这个值，可自行修改<br /> !undef为数据中没有这么多层的时候用来填补的缺测值<br /> !startname和endname分别为数据开始和结束的时间，用来动态生成文件名，可自行修改<br /> !程序执行后会分别生成theta和v的两个数据文件，每个文件中都有N个时次，注意该程序每次最多只能运行01时到23时，如果要更多时次，请自行修改程序中文件名的动态生成方法。<br /> integer,parameter::zcount=100,undef=32766,startname=2012081201,endname=2012081206<br /> real var(5)<br /> character cname*21<br /> integer::status=0,iname<br /> <br /> open(2,file='20120812.theta.grd',status='replace',form='binary')<br /> open(3,file='20120812.v.grd',status='replace',form='binary')<br /> <br /> do iname=startname,endname<br /> call getCname(iname,cname)<br /> open(1,file='file\'//cname,status='old')<br /> do i=1,zcount<br />   read(1,*,iostat=status)var<br />   if(status/=0)then<br />     write(2)undef<br />     write(3)undef<br />   else<br />     write(2)var(4)<br />     write(3)var(5)<br />   endif<br /> enddo<br /> close(1)<br /> enddo<br /> close(2)<br /> close(3)<br /> end<br /> <br /> subroutine getCname(iname,cname)<br /> integer::iname<br /> character*21 cname<br /> cname=''<br /> write(cname,'(i10)')iname<br /> cname=trim(cname)//'calwswd.dat'<br /> return<br /> endsubroutine<br /> <br />

生成数据文件之后增加两个ctl描述文件：

1. #v.ctl
   
2. dset ^20120812.v.grd
   
3. title outputv
   
4. undef 32766
   
5. XDEF  1 LINEAR 1 1
   
6. YDEF  1 LINEAR 1 1
   
7. ZDEF 100 linear 0 0.02
   
8. TDEF  6 LINEAR 00Z01AUG2012 1hr
   
9. vars 1
   
10. v 100 99 v data
    
11. endvars
    

#v.ctl<br /> dset ^20120812.v.grd<br /> title outputv<br /> undef 32766<br /> XDEF  1 LINEAR 1 1<br /> YDEF  1 LINEAR 1 1<br /> ZDEF 100 linear 0 0.02<br /> TDEF  6 LINEAR 00Z01AUG2012 1hr<br /> vars 1<br /> v 100 99 v data<br /> endvars<br />

1. #theta.ctl
   
2. dset ^20120812.theta.grd
   
3. title outputtheta
   
4. undef 32766
   
5. XDEF  1 LINEAR 1 1
   
6. YDEF  1 LINEAR 1 1
   
7. ZDEF 100 linear 0 0.02
   
8. TDEF  6 LINEAR 00Z01AUG2012 1hr
   
9. vars 1
   
10. theta 100 99 theta data
    
11. endvars
    

#theta.ctl<br /> dset ^20120812.theta.grd<br /> title outputtheta<br /> undef 32766<br /> XDEF  1 LINEAR 1 1<br /> YDEF  1 LINEAR 1 1<br /> ZDEF 100 linear 0 0.02<br /> TDEF  6 LINEAR 00Z01AUG2012 1hr<br /> vars 1<br /> theta 100 99 theta data<br /> endvars<br />

最后编写脚本作图：

1. 'reinit'
   
2. 'open C:\Users\fangbao\Desktop\test\theta.ctl'
   
3. 'open C:\Users\fangbao\Desktop\test\v.ctl'
   
4. 'set gxout barb'
   
5. 'set t 1 6'
   
6. 'set z 1 100'
   
7. 'd -sin(theta*3.14159/180)*v.2*2.5;-cos(theta*3.14159/180)*v.2*2.5;v.2'
   
8. 'printim C:\Users\fangbao\Desktop\test\test.png'
   
9. ;
   

'reinit'<br /> 'open C:\Users\fangbao\Desktop\test\theta.ctl'<br /> 'open C:\Users\fangbao\Desktop\test\v.ctl'<br /> 'set gxout barb'<br /> 'set t 1 6'<br /> 'set z 1 100'<br /> 'd -sin(theta*3.14159/180)*v.2*2.5;-cos(theta*3.14159/180)*v.2*2.5;v.2'<br /> 'printim C:\Users\fangbao\Desktop\test\test.png'<br /> ;<br />

有一个要注意的地方，GrADS中时间是00-23，注意和你的数据进行对应，由于仅做一个示例，我没有特别注意这一点
从你的数据可以看出，开始的数据特别大，绘制的图形也体现出来了（图形未做美化）：


整个使用的文件打包：
wnd.rar (29.92 KB, 下载次数: 20)

2012-10-15 13:32 上传

点击文件名下载附件
欢迎楼主多多介绍宝岛的同行加入~