- 积分
- 7124
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2017-11-26
- 最后登录
- 1970-1-1
|

楼主 |
发表于 2020-6-18 10:47:12
|
显示全部楼层
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; Plot data on a cross section
; This script will plot data at a set angle through a specified point
; This script adds lon/lat info along X-axis
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
dir = "/home/szw/wrfout/"
filename = "wrfout_d01_2020-06-04_12:00:00"
a = addfile(dir + filename + ".nc","r")
; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
type = "png"
;wks = gsn_open_wks(type,"plt_UW_mdBz_08.ncl")
;wks = gsn_open_wks(type,"plt_UW_mdBz_08.ab.2.ncl")
wks = gsn_open_wks(type,"plt_UW_mdBz_08.ab3.ncl")
;gsn_define_colormap(wks,"WhViBlGrYeOrReWh") ; Overwrite the standard color map
;gsn_define_colormap(wks, "Radar2") ;加载自己定义的色标
gsn_define_colormap(wks, "Radar.CINRAD") ;加载自己定义的色标
; Set some basic resources
res = True
res@MainTitle = "REAL-TIME WRF"
res@Footer = False
pltres = True
mpres = True
mpres@mpGeophysicalLineColor = "Black"
mpres@mpNationalLineColor = "Black"
mpres@mpUSStateLineColor = "Black"
mpres@mpGridLineColor = "Black"
mpres@mpLimbLineColor = "Black"
mpres@mpPerimLineColor = "Black"
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpOutlineSpecifiers = (/"China:states","Taiwan","Zhejiang"/)
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpNationalLineThicknessF = 2
mpres@mpNationalLineColor = "black"
mpres@mpGeophysicalLineColor = "black"
mpres@mpProvincialLineColor = "black"
mpres@mpProvincialLineThicknessF = 2
mpres@mpFillOn = True
mpres@mpFillAreaSpecifiers = (/"water", "land" /)
mpres@mpSpecifiedFillColors = (/"skyblue1","white"/)
mpres@mpMaskAreaSpecifiers = (/"China:states","Taiwan","Zhejiang"/)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTime = True
FirstTimeMap = True
times = wrf_user_getvar(a,"times",-1) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
dNames = getfilevardims(a,"P") ; NCL 6.4.0 and earlier
nd = dimsizes(mdims) ; 查看mdims为几围数组
print(dNames + " : " + mdims) ; 查看各个变量维数
xlat = wrf_user_getvar(a, "XLAT",0)
xlon = wrf_user_getvar(a, "XLONG",0)
ter = wrf_user_getvar(a, "HGT",0)
u10 = wrf_user_getvar(a,"U10",0) ; u at 10 m, mass point
v10 = wrf_user_getvar(a,"V10",0) ; v at 10 m, mass point
u10@units = "m/s"
v10@units = "m/s"
printVarSummary(xlat)
;---------------------------------------------------------------
do it = 0,ntimes-1 ; TIME LOOP
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
z = wrf_user_getvar(a, "z",it) ; grid point height
dbz = wrf_user_getvar(a,"dbz",it) ;dBz
mdbz = wrf_user_getvar(a,"mdbz",it) ;dBz
tc = wrf_user_getvar(a, "tc", it)
p = wrf_user_getvar(a, "pressure",it) ;
rh = wrf_user_getvar(a,"rh",it) ; relative humidity;
printVarSummary(tc)
u = wrf_user_getvar(a,"ua",it)
v = wrf_user_getvar(a,"va",it)
w = wrf_user_getvar(a,"wa",it)
if ( FirstTime ) then ; get height info for labels
zmin = 0.
zmax = 3. ; We are only interested in the first 6km
nz = floattoint(zmax + 1)
end if
;---------------------------------------------------------------
do ip = 1, 1,1 ; we are doing 3 plots;我们这里只绘制2张图
; all with the pivot point (plane) in the center of the domain
; at angles 0, 45 and 90
;
; |
; angle=0 is | ;垂直剖面直线AB的形状
; |
;
;线段AB的起点和终点经纬度设置
;08.ab UTC
startlon = 116
startlat = 36
endlon = 118
endlat = 37
;startlat = 29.35
;endlat = 28.6
PI=3.1415926
;--------------------------------------------------------------------------
;经纬度转为网格点后对应的最小最大值
;--------------------------------------------------------------------------
; Xstart = xy(0,0) ;经度最小
; Xend = xy(0,1) ;经度最大
; Ystart = xy(1,0) ;纬度最小
; Yend = xy(1,1) ;纬度最大
;找到dbz最大值的点
xy = wrf_user_ll_to_ij(a,(/startlon,endlon/),(/startlat,endlat/),res) ;AB经纬度转换为坐标点
latlon = wrf_user_ij_to_ll(a,xy(0,:),xy(1,:),res)
print("----")
print("Requested min/max xlon = " + startlon + "/" + endlon)
print("Actual min/max xlon = " + xlon(xy(1,0)-1,xy(0,0)-1) + "/" + \
xlon(xy(1,1)-1,xy(0,1)-1))
print("Calculated min/max xlon = " + latlon(0,0) + "/" + latlon(0,1))
print("----")
print("Requested min/max xlat = " + startlat + "/" + endlat)
print("Actual min/max xlat = " + xlat(xy(1,0)-1,xy(0,0)-1) + "/" + \
xlat(xy(1,1)-1,xy(0,1)-1))
print("Calculated min/max xlat = " + latlon(1,0) + "/" + latlon(1,1))
print("----")
xy = xy-1 ;转换为NCL下标
print("经度格点0 = "+ xy(0,0)) ;经度
print("经度格点1 = "+ xy(0,1)) ;经度
print("----")
print("纬度格点0 = "+ xy(1,0)) ;纬度
print("纬度格点1 = "+ xy(1,1)) ;纬度
print("----")
;--------------------------------------------------------------------------
;本次经纬度转为网格点后对应的最小最大值
;--------------------------------------------------------------------------
; startlon = xy(0,0) ;经度最小
; startlat = xy(1,0) ;纬度最大
; endlon = xy(0,1) ;经度最大
; endlat = xy(1,1) ;纬度最小
;plot AB垂直剖面的位置
plane = new(4,float) ; a new array
;plane = (/ xy(0,0), xy(1,1)/) ;A(x,y)点坐标
;plane = (/ xy(0,1), xy(1,0)/) ;B(x,y)点坐标
plane = (/ xy(0,0), xy(1,0), xy(0,1), xy(1,1) /) ;A(x,y)点到B(x,y)坐标
;alfa = atan2((endlat-startlat),(endlon-startlon)) ;计算线段AB与y轴方向的夹角
alfa = atan2((endlon-startlon),(endlat-startlat)) ;计算线段AB与y轴方向的夹角
angle = 90-alfa*180/PI
;uv = u*cos(alfa)+v*sin(alfa) ;由于atan2函数正北为0度,因此线段AB与x轴方向的夹角angle和atan2函数出来的alfa之间要换算过
uv = u*sin(alfa)+v*cos(alfa) ;由于atan2函数正北为0度,因此线段AB与x轴方向的夹角angle和atan2函数出来的alfa之间要换算过
print("angle of Line AB with X in plot is " +angle+" Degree")
;opts = False ; start and end point not supplied
opts = True ; 给定线段AB的话,这个要设置为True才行!!!
;--------------------------------------------------------------------------
;由于atan2函数正北为0度,因此线段AB的夹角angle和atan2函数出来的alfa之间要换算过
;--------------------------------------------------------------------------
if(ip .eq. 1) then
;angle = 90.
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) ;This function interpolates 2D model data onto a specified line.
X_desc = "longitude"
end if
if(ip .eq. 2) then
;angle = 0.
X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
X_desc = "latitude"
end if
if(ip .eq. 3) then
;angle = 45.
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
X_desc = "longitude"
end if
;把数据插值到各个高度上
dbz_plane = wrf_user_intrp3d(dbz,z,"v",plane,angle,opts)
u_plane = wrf_user_intrp3d(uv,z,"v",plane,angle,opts)
w_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts)
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
mdbz_plane = wrf_user_intrp2d(mdbz,plane,angle,opts)
p_plane = wrf_user_intrp3d( p,z,"v",plane,angle,opts)
rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts)
printVarSummary(tc_plane)
printVarSummary(rh_plane)
printVarSummary(p_plane)
print("Max terrain height in plot is " + max(ter_plane)+"m")
;*计算假相当位温
pressure_levels = (/ 1000.,950.,900.,850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,\
350.,300.,250.,200.,150.,100./) ; pressure levels to plot
nlevels = dimsizes(pressure_levels) ; number of pressure levels
tk = tc_plane+273.166 ; W8 ]: G. [% v
printVarSummary(tk)
printVarSummary(pressure_levels)
;按指定匹配维扩展到与指定变量同等大小
;假定数组q的维数 为 nt x ny x nx x nz ,而数组dz 为一维数组,长度为nz
;按 q 的第3维dz扩展 到 q 的大小
;Assume ntim1=100 and ntim2=1, and that T1 is a 4D array with dimensions ntim1 x klev x nlat x mlon and dp is an array with dimensions ntim2 x klev. Note that ntim2 is a degenerate dimension.
;In NCL V6.3.0 and later, if you want to conform dp to the same size of T1:
;dpC = conform_dims(dimsizes(T), dp, (/0,1/)) ; => (ntim1,klev,nlat,mlon)
P = p_plane
a1 = where(tk .gt. 263.0, 0.622*6.11*exp(17.26*(tk-273.16)/(tk-35.86)), \
0.622*6.11*exp(21.87*(tk-273.16)/(tk-7.66))) ;
b1= where(tk .gt. 263.0, P-0.278*exp(17.26*(tk-273.16)/(tk-35.86)),\
P-0.278*exp(21.87*(tk-273.16)/(tk-7.66))) ;
qs1=a1/b1
q1=qs1*rh_plane
e1=P*q1/100./(0.62197+q1/100.0)
tk1=55.0+2840.0/(3.5*log(tk)-log(e1)-4.805)
pot1=tk*(1000/P)^(0.2854*(1.0-0.28*q1/100.0))
ept1=pot1*exp(((3376./tk1)-2.54)*q1/100.0*(1.0+0.81*q1/100.0)) ;
ept1@description = "0se"
ept1@units = "K"
copy_VarCoords(tc_plane,ept1) ; assign coordinates
printVarSummary(ept1)
; Find the index where 6km is - only need to do this once
if ( FirstTime ) then
zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
b = ind(zz(:,0) .gt. zmax*1000. )
zmax_pos = b(0) - 1
if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then
zspan = b(0) - 1
else
zspan = b(0)
end if
delete(zz)
delete(b)
FirstTime = False
end if
; X-axis lables
dimsX = dimsizes(X_plane)
xmin = X_plane(0)
xmax = X_plane(dimsX(0)-1)
xspan = dimsX(0)-1
nx = floattoint( (xmax-xmin)/2 + 1)
;nx = 10
print("xmin and xmax in plot is " + xmin+" and "+xmax)
print("xspan and nx in plot is " + xspan+" and "+nx)
;---------------------------------------------------------------
; Options for XY Plots
opts_xy = res
opts_xy@tiXAxisString = X_desc
opts_xy@tiYAxisString = "Height (km)"
opts_xy@cnMissingValPerimOn = True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
opts_xy@tmXTOn = False
opts_xy@tmYROn = False
opts_xy@tmXBMode = "Explicit"
;;-------------------想要绘制的数据的经纬度范围-------------------------------------------------
minlat = startlat
maxlat = endlat
minlon = startlon
maxlon = endlon
loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),res) ;查找目标经纬度对应的格点
loc = loc-1 ;转换为NCL下标
; Xstart = loc(0,0) ;经度最小
; Xend = loc(0,1) ;经度最大
; Ystart = loc(1,0) ;纬度最小
; Yend = loc(1,1) ;纬度最大
; X-axis lables
xmin = minlon
xmax = maxlon
nx = 10
;;-----------------------------------------------------------------------------------------------
opts_xy@tmXBMode = "Explicit"
opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks
opts_xy@tmXBLabels = sprintf("%.2f",fspan(xmin,xmax,nx)) ; Create labels
;opts_xy@tmXBLabels = (/"29.60N,119.20E","29.40N,119.40E","29.20N,119.60E","29.00N,119.80E","28.80N,120.00E"/) ; Create labels
opts_xy@tmXBLabelFontHeightF = 0.010
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy@tiXAxisFontHeightF = 0.020
opts_xy@tiYAxisFontHeightF = 0.020
opts_xy@tmXBMajorLengthF = 0.02
opts_xy@tmYLMajorLengthF = 0.02
opts_xy@tmYLLabelFontHeightF = 0.015
; Plotting options for dBz
opts_dbz = opts_xy
opts_dbz@cnFillOn = True
;opts_dbz@ContourParameters = (/ 5., 65., 5./)
opts_dbz@cnLevelSelectionMode= "ManualLevels"
opts_dbz@cnMinLevelValF = 0.
opts_dbz@cnMaxLevelValF = 65.
opts_dbz@cnLevelSpacingF = 5.0 ;色标间隔
; Plotting options for mdBz
opts_mdbz = opts_xy
opts_mdbz@cnFillOn = True
;opts_mdbz@ContourParameters = (/ 5., 65., 5./)
opts_mdbz@cnLevelSelectionMode= "ManualLevels"
opts_mdbz@cnMinLevelValF = 0.
opts_mdbz@cnMaxLevelValF = 65.
opts_mdbz@cnLevelSpacingF = 5.0 ;色标间隔
; Plotting options for terrain
opts_ter = opts_xy
opts_ter@gsnYRefLine = 0.0
opts_ter@gsnAboveYRefLineColor = "black"
opts_ter@gsnDraw = False
opts_ter@gsnFrame = False
opts_ter@trYMaxF = zmax*1000 ;地形数据为m,因此这里要放大倍数为km
vcres = opts_xy
;vcres@vcGlyphStyle = "WindBarb"
vcres@vcGlyphStyle = "CurlyVector";"LineArrow";
vcres@vcPositionMode = "ArrowTail"
vcres@vcRefAnnoOn = True
vcres@vcMinDistanceF = 0.03
vcres@vcRefMagnitudeF = 0.5 ;单位长度
vcres@vcLineArrowThicknessF = 3.0
vcres@vcMapDirection = False
vcres@vcFillArrowsOn = True
vcres@vcRefLengthF = 0.015 ;矢量长度
vcres@vcFillArrowHeadXF = 0.2
vcres@vcFillArrowHeadYF = 0.2
vcres@vcFillArrowHeadInteriorXF = 0.5
;;u ,v for streamline
uvres = opts_xy ; plot mods desired
uvres@tiMainString = "Streamlines" ; title
uvres@stArrowLengthF = 0.004 ; size of the arrows.
uvres@stMinArrowSpacingF = 0.004 ; arrow spacing.
uvres@stArrowStride = 1 ; arrows start every third
uvres@stArrowLengthF = 0.008 ; default is dynamic
uvres@stLengthCheckCount = 15 ; default is 35
uvres@stLineStartStride = 3.5 ; default is 2 流线的密度
uvres@stMinArrowSpacingF = 0.035 ; default is 0.0
uvres@stStepSizeF = 0.001 ; default is dynamic
uvres@stLineThicknessF = 2.0 ; changes the line thickness
;;plot for tc温度
opts_tc = opts_xy
opts_tc@cnFillOn = False
opts_tc@cnLineColor="Black"
opts_tc@cnLineThicknessF = 2.0
;opts_tc@cnLineLabelFontColor="Black"
opts_tc@cnInfoLabelOrthogonalPosF=0.00
opts_tc@ContourParameters=(/ 2.5 /)
opts_tc@gsnContourPosLineDashPattern=1
opts_tc@cnLineDashPattern=1
;;plot for 0se相当位温
opts_0se = opts_xy
opts_0se@cnFillOn = False
opts_0se@cnLineColor="Blue"
opts_0se@cnLineThicknessF = 4.0
opts_tc@cnLineLabelFontColor="Blue"
opts_0se@cnInfoLabelOrthogonalPosF=0.00
opts_0se@ContourParameters=(/ 4.0 /)
opts_0se@gsnContourPosLineDashPattern=0
opts_0se@cnLineDashPattern=0
; Get the contour info for the rh and temp
contour_dbz = wrf_contour(a,wks,dbz_plane(0:zmax_pos,:),opts_dbz) ;;loc(0,0):loc(0,1)为想要绘制的经度范围
;contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)
;vector = gsn_streamline(wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:)*10.,uvres)
vector = gsn_vector(wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:)*10.,vcres)
contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc) ;;loc(0,0):loc(0,1)为想要绘制的经度范围
contour_0se = wrf_contour(a,wks,ept1(0:zmax_pos,:),opts_0se) ;;loc(0,0):loc(0,1)为想要绘制的经度范围
; contour_dbz = wrf_contour(a,wks,dbz_plane(0:zmax_pos,loc(0,0):loc(0,1)),opts_dbz) ;;loc(0,0):loc(0,1)为想要绘制的经度范围
contour_ter = gsn_csm_xy(wks,X_plane,ter_plane(loc(0,0):loc(0,1)),opts_ter)
; vector = gsn_streamline(wks,u_plane(0:zmax_pos,loc(0,0):loc(0,1)),w_plane(0:zmax_pos,loc(0,0):loc(0,1))*10.,uvres)
; MAKE PLOTS
if (FirstTimeMap) then
;lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
;lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
pltres = True
pltres@FramePlot = False
;mpres = True
optsM = res
optsM@NoHeaderFooter = True
optsM@cnFillOn = True
optsM@lbTitleOn = True;False
optsM@cnLevelSelectionMode= "ManualLevels"
optsM@cnMinLevelValF = 0.
optsM@cnMaxLevelValF = 65.
optsM@cnLevelSpacingF = 5.0 ;色标间隔
optsM@pmLabelBarOrthogonalPosF = -0.12 ;!!!!!!!!!!!!!色标与图的距离 !!!!!!!!!!!!!
;设置mdbz绘图区域
minlat=36
minlon=116
maxlat=37
maxlon=118
loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),res) ;查找目标经纬度对应的格点
contour = wrf_contour(a,wks,mdbz(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),optsM)
; Plotting options for Wind Vectors
opts = res
opts@FieldTitle = "Wind 10 m" ; overwrite Field Title
opts@NumVectors = 47 ; density of wind barbs
opts@vcWindBarbScaleFactorF = 2.5 ; 风杆符合国内习惯(4m/s一个长杆)
opts@vcRefMagnitudeF = 8. ; make vectors larger
opts@vcRefLengthF = 0.040 ; ref vec length
opts@vcGlyphStyle = "LineArrow" ; select wind barbs
opts@vcMinDistanceF = 0.025 ; thin out windbarbs
opts@vcLineArrowThicknessF =1.5
opts@vcRefAnnoOn = True ;风速单位标注
;;绘制水平风场
vector_uv = wrf_vector(a,wks,u10(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),v10(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),opts)
delete(opts)
;设置绘图区域
mpres@mpLeftCornerLatF = minlat
mpres@mpRightCornerLatF = maxlat
mpres@mpLeftCornerLonF = minlon
mpres@mpRightCornerLonF = maxlon
;plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
plot = wrf_map_overlays(a,wks,(/contour,vector_uv/),pltres,mpres)
;绘制垂直剖面的直线AB
lnres = True
lnres@gsLineThicknessF = 6.0
lnres@gsLineColor = "Black"
; do ii = 0,dimsX(0)-2
; gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres)
; end do
gsn_polyline(wks,plot,(/startlon,endlon/),(/startlat,endlat/),lnres)
frame(wks)
;delete(lon_plane)
;delete(lat_plane)
pltres@FramePlot = True
end if
plot = wrf_overlays(a,wks,(/contour_dbz,contour_tc,vector,contour_ter/),mpres) ; plot x-section
;plot = wrf_overlays(a,wks,(/contour_dbz,vector,contour_ter,contour_0se/),mpres) ; plot x-section
; Delete options and fields, so we don't have carry over
delete(opts_xy)
delete(opts_dbz)
delete(vcres)
delete(X_plane)
delete(ter_plane)
delete(mdbz_plane)
end do ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTimeMap = False
end do ; END OF TIME LOOP
end |
|