- 积分
- 46
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2022-7-20
- 最后登录
- 1970-1-1
|

楼主 |
发表于 2023-11-20 19:48:33
|
显示全部楼层
这是代码
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
a = addfile("/home/ychz/wrfout_d03_2022-07-03_00:00:00","r")
wks = gsn_open_wks("pdf","tem")
;---Set common resources for all plots
res = True
res@NoHeaderFooter = True ; Switch headers and footers off
res@pmLabelBarOrthogonalPosF = -0.1
res@lbTitleOn = False
pltres = True
pltres@PanelPlot = True ; Indicate these plots are to be paneled.
map_res = True
;res@gsnFrame = False
;res@gsnDraw = False
;res@gsnLeftString = ""
;res@gsnRightString = ""
; What times and how many time steps are in the data set?
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
ntimes = dimsizes(times) ; number of times in the file
plots = new ( 2, graphic )
do it = 0,ntimes-1,6 ; TIME LOOP
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
;---Read several WRF variables
tc = wrf_user_getvar(a,"tc",it) ; 3D temperature
u = wrf_user_getvar(a,"ua",it) ; 3D U at mass points
v = wrf_user_getvar(a,"va",it) ; 3D V at mass points
lat2d = wrf_user_getvar(a,"XLAT",it) ; latitude
lon2d = wrf_user_getvar(a,"XLONG",it) ; longitude
;---Now get the lowest (bottommost) level
nl = 0
tc2 = tc(nl,:,:)
u10 = u(nl,:,:)
v10 = v(nl,:,:)
tf2=tc2
w10 = sqrt(u10*u10+v10*v10)
u10 = u10*1.94386 ; Convert wind into knots
v10 = v10*1.94386
;---Change the metadata
tf2@description = "2m Temperature"
tf2@units = "degC"
u10@units = "kph"
v10@units = "kph"
;---Necessary for contours to be overlaid correctly on WRF projection
res@tfDoNDCOverlay = True ; Tell NCL you are doing a native plot
; res@tfDoNDCOverlay = "NDCViewport" ; can use this in NCL V6.5.0 or later
;---Temperature filled contour plot
opts = res
opts@cnFillOn = True
opts@cnLevelSelectionMode = "ExplicitLevels"
opts@cnLinesOn = False
opts@cnLevels = ispan(20,34,1)
opts@lbLabelFontHeightF = 0.015
;tf_res@lbOrientation = "Vertical"
opts@pmLabelBarOrthogonalPosF = -0.005
contour_tf = wrf_contour(a,wks,tf2,opts)
delete(opts)
;---Wind vector plot
opts = res
opts@vcMinDistanceF = 0.03
opts@vcRefLengthF = 0.02
opts@vcMinFracLengthF = 0.2
opts@vcGlyphStyle = "WindBarb"
opts@vcRefAnnoOn = False
vector = wrf_vector(a,wks,u10,v10,opts)
delete(opts)
;---wind filled contour plot
opts = res
opts@cnFillOn = True
opts@cnLevelSelectionMode = "ExplicitLevels"
opts@cnLinesOn = False
opts@cnLevels = ispan(0,4,1)
opts@lbLabelFontHeightF = 0.015
;w10_res@lbOrientation = "Vertical"
opts@pmLabelBarOrthogonalPosF = -0.005
contour_w10 = wrf_contour(a,wks,w10,opts)
delete(opts)
;---Map plot
map_res = True
map_res@gsnFrame = False
map_res@gsnDraw = False
map_res@tiMainString = "10mwind and 2mtem" ;name
map_res@gsnLeftString = tf2@description + " (" + tf2@units + ")~C~" + \
"10m Wind (" + u10@units + ")"
map_res@gsnLeftStringFontHeightF = 0.01
opts = True
opts@cnFillOn = True
opts@sfXArray = lon2d
opts@sfYArray = lat2d
pltres@LatLonOverlay = True
map_res@mpLeftCornerLatF = 30.3 这里不知道怎么改!!!
map_res@mpRightCornerLatF = 31.1
map_res@mpLeftCornerLonF = 103.6
map_res@mpRightCornerLonF = 104.5
;---Set map resources based on projection on WRF output file
map_res = wrf_map_resources(a,map_res)
map = gsn_csm_map(wks,map_res)
;---Overlay plots on map and draw.
;overlay(map,contour_tf)
;overlay(map,vector)
plots(0) = wrf_map_overlays(a,wks,(/contour_tf/),pltres,map_res)
plots(1) = wrf_map_overlays(a,wks,(/contour_w10,vector/),pltres,map_res)
sichuan_shp_name = "/home/ychz/map/sichuan/sichuan.shp"
chengdu_shp_name = "/home/ychz/map/sichuan/chengdu.shp"
lnres = True
lnres@gsLineColor = "black"
lnres@gsLineThicknessF = 1.5
sichuan_id = gsn_add_shapefile_polylines(wks,(/plots/),sichuan_shp_name,lnres)
chengdu_id = gsn_add_shapefile_polylines(wks,(/plots/),chengdu_shp_name,lnres)
pnlres = True
pnlres@txString = "PLOTS for : " + times(it)
pnlres@gsnPanelYWhiteSpacePercent = 3 ; Add white space b/w plots.
pnlres@gsnPanelScalePlotIndex = 1
gsn_panel(wks,(/plots/),(/2,1/),pnlres)
; draw(map) ; This will draw all overlaid plots and the map
frame(wks)
end do ; END OF TIME LOOP
end
|
|