激光测风雷达数据出图

wsx

begin
a = asciiread("CDL_3D10K_lidar3D10K-099_DBS5_PitchAngle71.38_Resolution015_StartIndex005_Ave10min_20250411 000010.csv",-1, "string")
yymmdd = new(dimsizes(a)-2,string)
hh     = new(dimsizes(a)-2,float)
minute = new(dimsizes(a)-2,float)
Temperature = new(dimsizes(a)-2,float)
Humidity = new(dimsizes(a)-2,float)
Pressure = new(dimsizes(a)-2,float)

;=======================================高度数据===========================================
level = new(298,float)
do j=0,297
    level(j) = stringtofloat(str_get_field(str_get_field(a(1), 5+j*9, ","),1,"m")) ;9:风速
end do
;======================================风数据读取============================================

wind_sped    = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_dir     = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_MeanSNR = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_DataObtainRate = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_StdDev  = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_MaxWindSpeed = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_MinWindSpeed = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_ZWind = new((/dimsizes(a)-2,298/),float);;;(time,level)
wind_ZWindStdDev = new((/dimsizes(a)-2,298/),float);;;(time,level)

do i = 0,dimsizes(a)-3
    a1 = a(i+2)
    a2 = str_get_field(a1, 1, ",")
    ;时间============================
    yymmdd(i) = str_get_field(a2, 1, " ")
    hh(i)     = stringtofloat(str_get_field(str_get_field(a2, 2 ," "),1,":"))
    minute(i) = stringtofloat(str_get_field(str_get_field(a2, 2 ," "),2,":"))
    ;变量============================
    Temperature(i) = stringtofloat(str_get_field(a1, 2, ","))
    Humidity(i) = stringtofloat(str_get_field(a1, 3, ","))
    Pressure(i) = stringtofloat(str_get_field(a1, 4, ","))
    ;风变量===========================
    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 5+j*9, ",")
    end do
    wind_sped(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 6+j*9, ",")
    end do
    wind_dir(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 7+j*9, ",")
    end do
    wind_MeanSNR(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 8+j*9, ",")
    end do
    wind_DataObtainRate(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 9+j*9, ",")
    end do
    wind_StdDev(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 10+j*9, ",")
    end do
    wind_MaxWindSpeed(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 11+j*9, ",")
    end do
    wind_MinWindSpeed(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 12+j*9, ",")
    end do
    wind_ZWind(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)

    wind_a4 = new(298,string)
    do j=0,297
        wind_a4(j) = str_get_field(a1, 13+j*9, ",")
    end do
    wind_ZWindStdDev(i,:) = stringtofloat(wind_a4)
    delete(wind_a4)
end do
;==========================测试-垂直速度=================================================
idx = ind(hh.ge.11.and.hh.le.23)

wind_ZWind_plot = transpose(wind_ZWind(idx,:)) ;绘图转秩
wind_ZWind_plot = where(abs(wind_ZWind_plot).gt.900, wind_ZWind_plot@_FillValue, wind_ZWind_plot)
time_plot = hh(idx) + minute(idx)/60.0

;==========================测试-信噪比=====================================================
wind_MeanSNR_plot = transpose(wind_MeanSNR(idx,:)) ;绘图转秩
wind_MeanSNR_plot = 10*log(wind_MeanSNR_plot)/log(10) ;量级处理

; 创建风速的等高线图
wks = gsn_open_wks("X11", "wind_profile")
res = True
res@cnFillOn = True
res@cnLevelSelectionMode = "ManualLevels"
;----------------------------------------
;res@cnMinLevelValF = -1
;res@cnMaxLevelValF = 1
;res@cnLevelSpacingF = 0.2
;----------------------------------------
res@cnMinLevelValF = 0
res@cnMaxLevelValF = 40
res@cnLevelSpacingF =2
res@cnLinesOn = False
res@sfXArray = time_plot         ; 指定X轴坐标;
res@sfYArray = level             ; 指定Y轴坐标
res@vpWidthF = 0.75       ; 图形宽度（占页面比例）
res@vpHeightF = 0.3
plot = gsn_csm_contour(wks, wind_MeanSNR_plot, res)
;========================================================================================



;==========================测试-风羽图=====================================================
wind_sped_plot = wind_sped(idx,:)
wind_dir_plot  = wind_dir(idx,:)

pi = 3.1415926
dir_rad = (270.0 - wind_dir_plot) * pi / 180.0
u = transpose(wind_sped_plot * cos(dir_rad))
v = transpose(wind_sped_plot * sin(dir_rad))
u = where(abs(u).gt.100, u@_FillValue, u)
v = where(abs(v).gt.100, v@_FillValue, v)
;-----------------------------------------------------------------------------------------
;-----------------------------------------------------------------------------------------
;===================================test=============================================================
vcres2 = True
vcres2@gsnDraw = False
vcres2@gsnFrame = False
vcres2@vcGlyphStyle = "WindBarb"
vcres2@vcWindBarbLineThicknessF = 1.0
vcres2@vcWindBarbColor = "black"
vcres2@vcRefLengthF = 0.04
vcres2@vcRefMagnitudeF = 10.0
vcres2@vcRefAnnoOn = False
vcres2@vcMinDistanceF = 0.02
vcres2@vcWindBarbScaleFactorF = 2.5
vcres2@vcPositionMode = "ArrowTail"
u!0 = "level"
u!1 = "time"
u&level = level
u&time  = time_plot

v!0 = "level"
v!1 = "time"
v&level = level
v&time  = time_plot
; ✅ 强制 Y 轴为线性（气压层通常是非等间距）
vcres2@gsnYAxisIrregular2Linear = True
vcres2@vpWidthF  = 0.75       ; 图形宽度（占页面比例）
vcres2@vpHeightF = 0.3


vcres2@vcMapDirection = False

plot2 = gsn_csm_vector(wks, u, v, vcres2)
draw(plot2)
frame(wks)

end