本篇主要介紹如何利用NCL軟體繪製颱風路徑圖，其中包括颱風路徑資料下載，NCL繪製颱風路徑，NCL繪製WRF模式模擬颱風路徑。
方法/步驟
颱風路徑資料下載
中國氣象局熱帶氣旋資料中心（
資料說明和資料下載如圖
開啟資料，格式如下，儲存為“typhoon.txt”
NCl指令碼如下
;********************************************************
; Load NCL scripts
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;********************************************************
begin
fiTY = "./typhoon.txt"
nrow = numAsciiRow(fiTY)
YYYYMMDDHH = new(nrow, "string")
lat = new(nrow, "float")
lon = new(nrow, "float")
vmax = new(nrow, "integer")
mslp = new(nrow, "integer")
data = asciiread(fiTY, -1, "string")
YYYYMMDDHH = str_get_field(data, 1," ")
lat = stringtofloat(str_get_field(data, 3," "))*0.1
lon = stringtofloat(str_get_field(data, 4, " "))*0.1
mslp = stringtoint(str_get_field(data, 5," " ))
vmax = stringtoint(str_get_field(data, 6, " "))
;print(vmax)
DateChar = stringtochar(YYYYMMDDHH)
MM = chartostring(DateChar(:,4:5))
DD = chartostring(DateChar(:,6:7))
HH = chartostring(DateChar(:,8:9))
HHi = mod(stringtoint(YYYYMMDDHH), 100)
DDi = mod(stringtoint(YYYYMMDDHH)/100, 100)
MMi = mod(stringtoint(YYYYMMDDHH)/10000, 100)
wks = gsn_open_wks("pdf", "suli_track")
gsn_define_colormap(wks,"rainbow")
res = True
[email protected] = False
[email protected] = False
[email protected] = 10
[email protected] = 40
[email protected] = 112
[email protected] = 160
[email protected] = 160
[email protected] = 0.025
[email protected] = 0.025
[email protected] = "True"
[email protected] = "MaskNotOcean"
[email protected] = 15
[email protected] = 5.0
[email protected] = True
[email protected] = "National"
[email protected] = "MediumRes"
[email protected] = "Earth..4"
[email protected] = "China:States"
plot = gsn_csm_map_ce(wks,res)
colours = (/3, 20, 60, 120, 180/)
resDot = True
resLine = True
dumDot= new(nrow, graphic)
dumLine = new(nrow, graphic)
do i = 0, nrow-2
xx = (/ lon(i), lon(i+1)/)
yy = (/ lat(i), lat(i+1)/)
if (vmax(i).le.17) then
[email protected] = colours(0)
end if
if (vmax(i) .ge. 17 .and. vmax(i).le.32) then
[email protected] = colours(1)
end if
if (vmax(i).ge.32 .and. vmax(i) .le. 42) then
[email protected] = colours(2)
end if
if (vmax(i).ge.42 .and. vmax(i) .lt. 51) then
[email protected] = colours(3)
end if
if (vmax(i).ge.51) then
[email protected] = colours(4)
end if
dumLine(i) = gsn_add_polyline(wks, plot, xx, yy, resLine)
end do
[email protected] = "black"
do i = 0, nrow-1
if (HH(i) .eq. "00") then
[email protected] = 0.02
dumDot(i) = gsn_add_polymarker(wks, plot, lon(i), lat(i), resDot)
else
[email protected] = 0.005
dumDot(i) = gsn_add_polymarker(wks, plot, lon(i), lat(i), resDot)
end if
end do
resLg = True
[email protected] = "MarkLines"
[email protected] = True
[email protected] = colours
[email protected] = 0.02
[email protected] = True
[email protected] = colours
[email protected] = 0.14
[email protected] = 0.13
[email protected] = "Background"
[email protected] = 0.3
[email protected] = 0.015
[email protected] = "Best Track"
lbid = gsn_create_legend(wks, 5, (/" 17m/s or less"," 17 to 32m/s"," 32 to 42m/s"," 42 to 51m/s"," 51 or more"/), resLg)
amres = True
[email protected] = 0.3
[email protected] = -0.3
dumLg = gsn_add_annotation(plot, lbid, amres)
dumDate = new(nrow,graphic)
resTx = True
[email protected] = 0.015
[email protected] = "black"
[email protected] = "BottomLeft"
do i = 0, nrow-1
if (HH(i) .ne. "00" ) then
continue
end if
dumDate(i) = gsn_add_text(wks,plot, MM(i)+DD(i), lon(i)+0.5, lat(i), resTx)
end do
draw(wks)
frame(wks)
end
如果用WRF輸出資料繪製颱風路徑，得到下圖
********************************************************
; Plot storm stracks from wrfout files.
;********************************************************
; ===========================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
; DATES
date = (/1512,1600,1612,1700,1712,1800,1812,1900/)
ndate = dimsizes(date)
sdate = sprinti("%4.0i",date)
; Experiment name (for legend)
EXP = (/"EXP_I"/) ; (/"EXP_I","EXP_II","EXP_III"/)
nexp = dimsizes(EXP)
; To get lat/lon info.
a = addfile("wrfout_d01_2003-07-15_00:00:00.nc","r")
lat2d = a->XLAT(0,:,:)
lon2d = a->XLONG(0,:,:)
dimll = dimsizes(lat2d)
nlat = dimll(0)
mlon = dimll(1)
; Sea Level Pressure
slp = wrf_user_getvar(a,"slp",0)
dims = dimsizes(slp)
; Array for track
time = new(ndate,string)
imin = new(ndate,integer)
jmin = new(ndate,integer)
smin = new(ndate,integer)
; =======
; ndate
; =======
fs = systemfunc("ls wrfout*00")
nfs= dimsizes(fs)
if(nfs .ne. ndate) then
print("Check input data:"+nfs+" .ne. "+ndate)
end if
do ifs=0,nfs-1
f = addfile(fs(ifs)+".nc","r")
time(ifs) = wrf_user_list_times(f)
; print(time(ifs))
slp2d = wrf_user_getvar(f,"slp",0)
; We need to convert 2-D array to 1-D array to find the minima.
slp1d = ndtooned(slp2d)
smin(ifs) = minind(slp1d)
; Convert the index for 1-D array back to the indeces for 2-D array.
minij = ind_resolve(ind(slp1d.eq.min(slp2d)),dims)
imin(ifs) = minij(0,0)
jmin(ifs) = minij(0,1)
; print(time(ifs)+" : "+min(slp2d)+" ("+imin(ifs)+","+jmin(ifs)+")")
end do
;
; Graphics section
wks=gsn_open_wks("ps","track") ; Open PS file.
gsn_define_colormap(wks,"BlGrYeOrReVi200") ; Change color map.
res = True
[email protected] = False ; Turn off draw.
[email protected] = False ; Turn off frame advance.
[email protected] = True ; Maximize plot in frame.
[email protected] = "Hurricane Isabel" ; Main title
WRF_map_c(a,res,0) ; Set up map resources
; (plot options)
plot = gsn_csm_map(wks,res) ; Create a map.
; Set up resources for polymarkers.
gsres = True
[email protected] = 16 ; filled dot
;[email protected] = 0.005 ; default - 0.007
cols = (/5,160,40/)
; Set up resources for polylines.
res_lines = True
[email protected] = 3. ; 3x as thick
dot = new(ndate,graphic) ; Make sure each gsn_add_polyxxx call
line = new(ndate,graphic) ; is assigned to a unique variable.
; Loop through each date and add polylines to the plot.
do i = 0,ndate-2
[email protected] = cols(0)
xx=(/lon2d(imin(i),jmin(i)),lon2d(imin(i+1),jmin(i+1))/)
yy=(/lat2d(imin(i),jmin(i)),lat2d(imin(i+1),jmin(i+1))/)
line(i) = gsn_add_polyline(wks,plot,xx,yy,res_lines)
end do
lon1d = ndtooned(lon2d)
lat1d = ndtooned(lat2d)
; Loop through each date and add polymarkers to the plot.
do i = 0,ndate-1
print("dot:"+lon1d(smin(i))+","+lat1d(smin(i)))
[email protected] = cols(0)
dot(i)=gsn_add_polymarker(wks,plot,lon1d(smin(i)),lat1d(smin(i)),gsres)
end do
; Date (Legend)
txres = True
[email protected] = 0.015
[email protected] = cols(0)
txid1 = new(ndate,graphic)
; Loop through each date and draw a text string on the plot.
do i = 0, ndate-1
[email protected] = "CenterRight"
ix = smin(i) - 4
print("Eye:"+ix)
if(i.eq.1) then
[email protected] = "CenterLeft"
ix = ix + 8
end if
txid1(i) = gsn_add_text(wks,plot,sdate(i),lon1d(ix),lat1d(ix),txres)
end do
; Add marker and text for legend. (Or you can just use "pmLegend" instead.)
[email protected] = "CenterLeft"
txid2 = new(nexp,graphic)
pmid2 = new(nexp,graphic)
do i = 0,nexp-1
[email protected] = cols(i)
[email protected] = cols(i)
ii = ((/129,119,109/)) ; ilat
jj = ((/110,110,110/)) ; jlon
ji = ii*mlon+jj ; col x row
pmid2(i) = gsn_add_polymarker(wks,plot,lon1d(ji(i)),lat1d(ji(i)),gsres)
txid2(i) = gsn_add_text(wks,plot,EXP(i),lon1d(ji(i)+5),lat1d(ji(i)),txres)
end do
draw(plot)
frame(wks)
end
