风廓线雷达数据出图

wsx

begin
time_band  = 30

a = asciiread("风廓线数据.txt",-1, "string")
data_uv_str = new((/time_band,54/),float) ;timeband=6
data_uv_dir = new((/time_band,54/),float) ;timeband=6
data_w_str  = new((/time_band,54/),float) ;timeband=6
data_w_CN2  = new((/time_band,54/),float) ;timeband=6
time_hour   = new((/time_band/),float) ;timeband=6
level_y     = new((/54/),float) ;timeband=6
do hour_idx = 0,time_band-1
;------------------数据读取------------------
i = 20+78*hour_idx ;每个时次54数据
yymmdd = str_get_field(a(i), 2, ":")
hour   = str_get_field(a(i+1), 2, ":")
level  = str_get_field(a(i+3:i+56), 2, " ")
uv_str = str_get_field(a(i+3:i+56), 3, " ")
uv_dir = str_get_field(a(i+3:i+56), 4, " ")
w_str  = str_get_field(a(i+3:i+56), 5, " ")
w_CN2  = str_get_field(a(i+3:i+56), 6, " ")
;-------------------------------------------
time_hour(hour_idx)     = stringtofloat(hour)
data_uv_str(hour_idx,:) = stringtofloat(uv_str)
data_uv_dir(hour_idx,:) = stringtofloat(uv_dir)
data_w_str(hour_idx,:)  = stringtofloat(w_str)
data_w_CN2(hour_idx,:)  = stringtofloat(w_CN2)
level_y                 = stringtofloat(level)
delete([/yymmdd,hour,level,uv_str,uv_dir,w_str,w_CN2,i/])
end do
;=========================绘图处理================
dir_rad = data_uv_dir * 3.14159 / 180.0
; 计算U和V分量（注意：风向是指风的来向）
; 气象上：U = -风速 * sin(风向), V = -风速 * cos(风向)
U = new((/time_band,54/),float)
V = new((/time_band,54/),float)
U = -data_uv_str * sin(dir_rad)
V = -data_uv_str * cos(dir_rad)
U = where(abs(U).gt.6000, U@_FillValue, U)
V = where(abs(V).gt.6000, V@_FillValue, V)
; 转置：将xy轴交换 (ntime, nlev) -> (nlev, ntime)
data_uv_str_T = transpose(data_uv_str)
U_T = transpose(U)
V_T = transpose(V)
;====================================================
; 绘制风羽图;已验证
wks = gsn_open_wks("X11", "风羽图")
vcres = True
vcres@gsnDraw = False
vcres@gsnFrame = False
vcres@vcGlyphStyle = "WindBarb"
vcres@vcRefLengthF = 0.04
vcres@vcMinDistanceF = 0.03
vcres@vcWindBarbColor = "black"
vcres@tmXBMode = "Explicit"
vcres@tmXBValues = ispan(0, time_band-1,10)
vcres@tmXBLabels = ispan(10,12,1)
vcres@tmYLMode = "Explicit"
vcres@tmYLValues = ispan(0, 53 ,1)
vcres@tmYLLabels = level_y
vector = gsn_csm_vector(wks, U_T, V_T, vcres)
draw(vector)
frame(wks)



end