求一个画WRFout台风路径的NCL脚本，不要官网的，要自己写的好用的！！谢谢！

vuchiyuan

利用WRF模拟台风，想用NCL画台风路径。
求一个画wrfout_d0*_2012_**_**_**:00:00台风路径的NCL脚本，不要官网的，要自己写的、好用的，最好带注释的！！谢谢！
下面这两个是官网给出的脚本及应用，改了几天了，一直出错，一方面赶时间，另外也实在没耐心了，因此来求一个。先谢啦～
官网链接如下，有耐心的童鞋可以试试。
http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/Examples/SPECIAL/wrf_Vortex.htm
http://www.ncl.ucar.edu/Applications/wrftrack.shtml

sun_shine_Xia

1. ;*********************************************
   
2. ; route.ncl
   
3. ;*********************************************
   
4. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
5. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   
6. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
   
7. load "/nuist/scratch/Xiayang/regcm_output/work/picture/cnmap/cnmap.ncl"
   
8. 
   
9. begin
   
10. 
    
11.   lines = new((/29,4/),"float")
    
12.   
    
13.   lines = asciiread("real_route.txt",(/29,4/),"float") ;以字符串形式读取参数文件入数组argu
    
14.    t   = lines(:,0)                                     ;从数组lines中获取时间
    
15.   lat  = lines(:,1)                                     ;从数组lines中获取经度值lon
    
16.   lon  = lines(:,2)                                     ;从数组lines中获取纬度值lat
    
17.   cslp = lines(:,3)
    
18.   times= flt2string(t)
    
19.   ;print(times)
    
20. 
    
21. ;***************** DATES ******************************
    
22.   ndate = dimsizes(t)
    
23.   
    
24.     f   = addfile("../wrfout_d02_2014-09-14_00:00:00.nc","r")
    
25.   lat2d = f->XLAT(0,:,:)
    
26.   lon2d = f->XLONG(0,:,:)
    
27.    times1 = wrf_user_getvar(f,"times",-1)  ; get all times in the file
    
28.   print(times1)
    
29.    slp  = wrf_user_getvar(f,"slp",-1)
    
30.    dims = dimsizes(slp(0,:,:))
    
31.   printVarSummary(slp)
    
32. 
    
33. ; Array for track
    
34.    time = new(79,string)
    
35.    imin = new(79,integer)
    
36.    jmin = new(79,integer)
    
37.    smin = new(79,integer)
    
38.    slpmin = new(79,float)
    
39. ;***************************
    
40. ;  ndate
    
41. ; ***************************
    
42. 
    
43.   do ifs=0,78
    
44.     slp2d = wrf_user_getvar(f,"slp",ifs)
    
45. 
    
46. ; We need to convert 2-D array to 1-D array to find the minima.
    
47.     slp1d     = ndtooned(slp2d)
    
48.     smin(ifs) = minind(slp1d)
    
49.     slpmin(ifs)= min(slp1d)
    
50. ; Convert the index for 1-D array back to the indeces for 2-D array.
    
51.     minij     = ind_resolve(ind(slp1d.eq.min(slp2d)),dims)
    
52.     imin(ifs) = minij(0,0)
    
53.     jmin(ifs) = minij(0,1)
    
54.   end do
    
55. 
    
56. ;***********************************************
    
57.   wks = gsn_open_wks("pdf","route2")               
    
58.      gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
    
59.   
    
60.   res                   = True
    
61.   res@gsnDraw           = False    ; will panel later
    
62.   res@gsnFrame          = False    ; will panel later
    
63. 
    
64.   res@mpLimitMode       = "LatLon"  
    
65.   res@mpMinLatF         = 12
    
66.   res@mpMaxLatF         = 28
    
67.   res@mpMinLonF         = 95
    
68.   res@mpMaxLonF         = 125
    
69.   res@mpCenterLonF      = 110
    
70.   res@gsnMajorLatSpacing= 4
    
71.   res@gsnMajorLonSpacing= 5
    
72.   res@gsnMinorLatSpacing= 2
    
73.   res@gsnMinorLonSpacing= 5
    
74.   res@tmXBLabelFontHeightF = 0.0175
    
75.   res@tmYLLabelFontHeightF = 0.0175
    
76.   
    
77.   
    
78.   ;res@gsnAddCyclic      = False
    
79. 
    
80.   res@mpFillOn                = False
    
81.   res@mpFillDrawOrder         = "PreDraw"
    
82.   res@mpOutlineOn             = True
    
83.   res@mpDataBaseVersion       = "MediumRes"
    
84.   res@mpDataSetName           = "Earth..4"
    
85.   ;res@mpAreaMaskingOn         = True
    
86.   ;res@mpMaskAreaSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
    
87.   res@mpLandFillColor         = -1
    
88.   res@mpInlandWaterFillColor  = -1
    
89.   res@mpOceanFillColor        = -1
    
90.   ;res@mpOutlineBoundarySets   = "NoBoundaries"
    
91. 
    
92.   cnres           = True
    
93.   cnres@china     = True       ;draw china map or not
    
94.   cnres@river     = False       ;draw changjiang&huanghe or not
    
95.   cnres@province  = True       ;draw province boundary or not
    
96.   cnres@nanhai    = False       ;draw nanhai or not
    
97.   cnres@diqu      = False       ; draw diqujie or not
    
98. 
    
99. 
    
100.   plot = gsn_csm_map_ce(wks,res)
     
101. 
     
102. 
     
103. ;***********************************************
     
104. ;            添加观测台风中心标示
     
105. ;***********************************************
     
106.    txres                 = True
     
107.    txres@txFont          = 137
     
108.    txres@txFontHeightF   = 0.025
     
109.    txres@txFontThicknessF= 4.0
     
110.    txres@txFontColor     = "black"
     
111.   do i = 10,28,4
     
112.     text=gsn_add_text(wks,plot,"p",lon(i),lat(i),txres)
     
113.   delete(text)
     
114.   end do
     
115.     text=gsn_add_text(wks,plot,"p",lon(28),lat(28),txres)
     
116.    
     
117. ;************************************************
     
118. ; add the Line
     
119. ;************************************************
     
120.   resp                  = True                      ; polyline mods desired
     
121.   resp@gsLineColor      = "black"                     ; color of lines
     
122.   resp@gsLineThicknessF = 5.0                       ; thickness of lines
     
123.   resp@gsLineLabelString= ""
     
124.   
     
125.     gsn_polyline(wks,plot,lon(10:28),lat(10:28),resp)     
     
126. 
     
127.   delete(resp)
     
128.   
     
129.   
     
130. ;************************************************
     
131. ; add the Time
     
132. ;************************************************  
     
133.    tres                 = True
     
134.    tres@txFont          = 25
     
135.    tres@txFontHeightF   = 0.0175
     
136.    tres@txFontThicknessF= 5
     
137.     do i = 10,28,4                 
     
138.         text=gsn_add_text(wks,plot,times(i),lon(i),lat(i)+1.8,tres)
     
139.     delete(text)
     
140.     end do
     
141.         text=gsn_add_text(wks,plot,times(28),lon(28)-1,lat(28)+1.8,tres)
     
142.     delete(text)
     
143.    
     
144. ;************************************************
     
145. ; add the (a)
     
146. ;************************************************  
     
147.    tres@txFontHeightF   = 0.0225
     
148.    text=gsn_add_text(wks,plot,"(a)",96,27,tres)
     
149. 
     
150. ;************************************************
     
151. ; add the simulaton line
     
152. ;************************************************   
     
153.   
     
154.    resp                  = True                      ; polyline mods desired
     
155.    resp@gsLineColor      = "red"                     ; color of lines
     
156.    resp@gsLineThicknessF = 5.0                       ; thickness of lines
     
157.    resp@gsLineLabelString= ""
     
158.    
     
159.    tres@txFontHeightF   = 0.02
     
160.    tres@txFontColor     = "red"
     
161.    txres@txFontColor    = "red"
     
162.      j=10
     
163.   do i = 24,77
     
164.      xx=(/lon2d(imin(i),jmin(i)),lon2d(imin(i+1),jmin(i+1))/)
     
165.      yy=(/lat2d(imin(i),jmin(i)),lat2d(imin(i+1),jmin(i+1))/)
     
166.        gsn_polyline(wks,plot,xx,yy,resp)
     
167.      
     
168.      if(.not.mod(i,12))then
     
169.        text1 = gsn_add_text(wks,plot,times(j),lon2d(imin(i),jmin(i)),lat2d(imin(i),jmin(i))-2,tres)
     
170.       
     
171.        text2 = gsn_add_text(wks,plot,"p",lon2d(imin(i),jmin(i)),lat2d(imin(i),jmin(i)),txres)   
     
172.        j=j+4
     
173.        delete(text1)
     
174.        delete(text2)     
     
175.      end if
     
176.   end do
     
177.        i = 78
     
178.        text1 = gsn_add_text(wks,plot,times(28),lon2d(imin(i),jmin(i))-1,lat2d(imin(i),jmin(i))-2,tres)
     
179.        text2 = gsn_add_text(wks,plot,"p",lon2d(imin(i),jmin(i)),lat2d(imin(i),jmin(i)),txres)
     
180. 
     
181.    print(slpmin)  
     
182.   
     
183.    write_table("simu_route.txt","w",[/times1(24),lat2d(imin(24),jmin(24)),lon2d(imin(24),jmin(24)),slp(24,imin(24),jmin(24))/],"%s %6.1f %6.1f %7.1f")
     
184.   do i=0,78,3
     
185.    write_table("simu_route.txt","a",[/times1(i),lat2d(imin(i),jmin(i)),lon2d(imin(i),jmin(i)),slp(i,imin(i),jmin(i))/],"%s %6.1f %6.1f %7.1f")
     
186.   end do
     
187.       
     
188. ;************************************************
     
189. ; add the China map
     
190. ;************************************************   
     
191.   chinamap = add_china_map(wks,plot,cnres)
     
192. 
     
193.   tres@txFontHeightF    = 0.0165
     
194.   txres@txFontHeightF   = 0.0172
     
195.   txres@txFontThicknessF= 3.0      
     
196.     lgtext1  = gsn_add_text(wks,plot,"p",99,16,txres)
     
197.     lgtext11 = gsn_add_text(wks,plot,"Simulation",102,16,tres)
     
198.    
     
199.   txres@txFontColor     = "black"
     
200.   tres@txFontColor      = "black"
     
201.     lgtext2 = gsn_add_text(wks,plot,"p",99,14.3,txres)
     
202.     lgtext22 = gsn_add_text(wks,plot,"Observation",102.2,14.3,tres)         
     
203.          
     
204.   draw(plot)
     
205.   frame(wks)
     
206. end

复制代码

vuchiyuan

sun_shine_Xia 发表于 2017-5-25 23:50
画出来就是上面这样的。红色的就是从模式读出来的台风中心路径

请问：
(1)第11行中，29,4是什么意思？
(2)第13行中，real_route.txt是什么文件？
(3)第34-38行中，79是什么意思？
谢谢！

vuchiyuan

sun_shine_Xia 发表于 2017-5-25 23:50
画出来就是上面这样的。红色的就是从模式读出来的台风中心路径

谢谢，这几个数基本上都猜出来了