Dear All,
I tried to mark the obs location using gsn_polymarker, and tried to
indicate a subdomain of interest using gsn_polyline. They work if I
don't add contour2 (using filled contour). However, if I add contour2
to the plot, then both the polymarker and the box don't show on the
plot. Please help me to find a solution. THX.
Yonghui
Following is the original code: in the plot I have two choice, one
without contour2, the other with contour2.
; Example script to produce an OLR plot for a WRF real-data run,
; with the ARW coordinate dynamics option.
; In this example we are also going to zoom into an area of interest
load "/opt/ncl_ncarg-5.1.1-pgi/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "/opt/ncl_ncarg-5.1.1-pgi/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "/opt/ncl_ncarg-5.1.1-pgi/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "/opt/ncl_ncarg-5.1.1-pgi/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "/opt/ncl_ncarg-5.1.1-pgi/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; model reference states
dir0 = "/raid2/yonghui/osse-exp/Result/fcst/True/"
datem = "10_18:00"
fname = "wrfout_d01_2008-02-" + datem + ":00.nc"
; obs network dimension
obsnorth = 34
obseast = 29
obssigma = 36
; gain information
Ibdy = 13 ;get rid of boundary
Lskip = 2 ;get Kalman gain every 2 grid
; dimension of each gain matrix 25*25
Ndim = 25
; obs location in model grid
Ix = 35 ;model grid in east
Iy = 77 ;model grid in north
Ilev = 7 ; vertical level
Ilev1 = Ilev - 1
; small domain centered by (Ix,Jy) with 2*Ibdy-1 as length for plot
x_start = Ix - Ibdy
x_end = Ix + Ibdy -2
y_start = Iy - Ibdy
y_end = Iy + Ibdy -2
; transform model grid coordinate of obs location in obs network grid
Ilon=(Ix-13)/2+1 ;obs network grid in east
Jlat=(Iy-13)/2+1 ;obs network grid in north
; record number
krec = ((Jlat-1)*obseast+Ilon-1)*obssigma+Ilev-1 ;since ncl start from
0 so -1
; gain data directory
pathd = "/raid2/yonghui/Gain/letkf-jump02/polynomial-dart-wholemean/"
; case name
dateh = "1018Z"
cases = "test-"+dateh+"-ens93-ndim8-nocenter-cut/"
; variable name
var="t"
; var="u"
; var="v"
; var="q"
vfile = "gain-"+var+"_org.bin"
; vfile = "gain-"+var+"_ret.bin"
; WRF input files needs to have a ".nc" appended, so just do it.
f0 = addfile(dir0 + fname,"r")
; read the gain
gaint = new ((/Ndim,Ndim/),"float",-9999.)
gaint = fbindirread(pathd+cases+vfile,krec,(/Ndim,Ndim/),"float")
; generate plots, but what kind do we prefer?
type = "pdf"
; type = "x11"
; type = "ps"
; type = "ncgm"
; output filename
wks = gsn_open_wks(type,"gain-"+dateh+"-Ix"+Ix+"-Iy"+Iy+"-lev"+Ilev)
; Change color map to something that has a grey scale
; gsn_define_colormap(wks,"wxpEnIR")
; Set some Basic Plot Information
res = True
; res@MainTitle = "OSSE" ; name of plot
res@InitTime = False ; do not plot initial time
res@Footer = False ; switch footers off
res@tiXAxisOn = False
res@tiYAxisOn = False
res@tiMainOn = False
res@tmXBLabelFontHeightF = 0.04 ; boost the axis value labels even more
res@tmYLLabelFontHeightF = 0.04
res@tmYRLabelFontHeightF = 0.04
; Set up two sets of plot and map resources
pltres = True
pltres@PanelPlot = True ; if want to put plots on 1 page
pltres@gsnDraw = False
pltres@gsnFrame = False
pltres@gsnMaximize = True
pltres@FramePlot = False ; donot advance frame automatically
pltres@lbLabelBarOn = False ; turn off label bar
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = wrf_user_list_times(f0) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
; map set
mpres = True
mpres@ZoomIn = True
mpres@Xstart = x_start
mpres@Ystart = y_start
mpres@Xend = x_end
mpres@Yend = y_end
mpres@tmXBLabelFontHeightF = 0.02 ; change lat lon font size
mpres@tmXBMajorLengthF = 0.02 ; change the tickmark length
mpres@mpGridAndLimbOn = True ; default is every 15 deg
mpres@mpGridSpacingF = 2 ; change to match labels
mpres@mpGeophysicalLineColor = "Gray"
mpres@mpNationalLineColor = "Gray"
mpres@mpUSStateLineColor = "Gray"
mpres@mpGridLineColor = "Black"
mpres@mpLimbLineColor = "Gray"
mpres@mpPerimLineColor = "Gray"
mpres@mpGeophysicalLineThicknessF = 1.5
mpres@mpGridLineThicknessF = 0.5
mpres@mpLimbLineThicknessF = 1.5
mpres@mpNationalLineThicknessF = 1.5
mpres@mpUSStateLineThicknessF = 1.5
mpres@gsnMaximize = True ; maximize image
t0 = wrf_user_getvar(f0,"tc",-1) ; get variable
u0 = wrf_user_getvar(f0,"ua",-1) ; get variable
v0 = wrf_user_getvar(f0,"va",-1) ; get variable
dimvar = dimsizes(t0)
print ("dimvar=" + dimvar)
print ("y_start=" + y_start)
print ("y_end=" + y_end)
print ("x_start=" + x_start)
print ("x_end=" + x_end)
Ilev1 = Ilev - 1 ;ncl start from 0 so -1
tzoom1 = t0(:,Ilev1,y_start:y_end,x_start:x_end) ; zoomed area
uzoom1 = u0(:,Ilev1,y_start:y_end,x_start:x_end) ; zoomed area
vzoom1 = v0(:,Ilev1,y_start:y_end,x_start:x_end) ; zoomed area
plot = new (1, graphic)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; loop over time
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do it = 0,ntimes-1 ; start at time 2
; contour plot of T
opts = res
opts@lbTitleOn = False ; remove field name from
label bar
opts@ContourParameters = (/-30, 30, 2/)
opts@lbLabelBarOn = False ; turn off label bar
opts@TimeLabel = times(it) ; set Valid time on the
plot
opts@cnLineColor = "black"
opts@cnInfoLabelOn = False ; no bottom info plotted
opts@pmLabelBarOrthogonalPosF = -0.1
opts@gsnContourLineThicknessesScale = 2.5
opts@cnLineLabelsOn = True
opts@gsnDraw = False ; Forces the plot to be drawn
opts@gsnFrame = False ; Frame advance
opts@mpGridAndLimbDrawOrder = "PreDraw"
opts@cnLabelDrawOrder = "PreDraw"
contour1 = wrf_contour(f0,wks,tzoom1(it,:,:),opts)
; replace t with the gain
dimzoom = dimsizes(tzoom1)
print ("dim of t0zoom=" + dimzoom)
tzoom1(it,:,:) = (/gaint(:,:)/)
optg = res
optg@cnInfoLabelOn = False
optg@UnitLabel = "C, shaded"
optg@FieldTitle = "Tdif" ; overwrite Field Title
; opts@gsnStringFontHeightF = 2
optg@lbLabelBarOn = True ; turn off label bar
optg@cnFillOn = True
optg@pmLabelBarOrthogonalPosF = -0.02
optg@ContourParameters = (/-1.2, 1.8, 0.3/)
optg@cnFillColors = (/"Blue3","Blue2","Blue1",\
"Cyan","SpringGreen","White","White","Yellow","Orange",\
"Red","Red2", "Red4"/)
optg@gsnContourLineThicknessesScale = 1.0
optg@gsnDraw = False ; Forces the plot to be drawn
optg@gsnFrame = False ; Frame advance
optg@lbTitleOn = False ; remove field name from
label bar
optg@TimeLabel = times(it) ; set Valid time on the
plot
optg@cnLineColor = "black"
optg@cnInfoLabelOn = False ; no bottom info plotted
optg@pmLabelBarOrthogonalPosF = -0.1
optg@gsnContourLineThicknessesScale = 2.5
optg@cnLineLabelsOn = True
optg@gsnDraw = False ; Forces the plot to be drawn
optg@gsnFrame = False ; Frame advance
contour2 = wrf_contour(f0,wks,tzoom1(it,:,:),optg)
; vector plot of (U,V)
resu = True
resu@TimeLabel = times(it) ; set Valid time on the
plot
resu@gsnMaximize = True ; Maximize plot in frame
resu@gsnSpreadColors = True ; span full colormap
resu@vcMonoLineArrowColor = True ; color arrows based on magnitude
resu@vcRefLengthF = 0.03313608
resu@vcMinFracLengthF = 0.3
;
; Setting the reference magnitude also affects the length
; of the arrows. In this case it is an inverse relationship.
;
resu@vcRefMagnitudeF = 50.0
;
; Line-drawn arrowheads are sized proportionally to the arrow length
; unless the resulting size would be outside the limits defined by
; the arrowhead minimum and maximum size resources. Setting the
; minimum and maximum sizes to the same value causes all the
;
resu@tiMainOn = False
resu@tiMainString = ""
resu@vcLineArrowHeadMinSizeF = 0.013
resu@vcLineArrowHeadMaxSizeF = 0.013
resu@gsnDraw = False ; don't draw
resu@gsnFrame = False ; don't advance frame
resu@vcExplicitLabelBarLabelsOn = False
resu@vcRefMagnitudeF = 15. ; make vectors larger
resu@vcRefLengthF = 0.050 ; reference vector length
resu@vcGlyphStyle = "CurlyVector" ; turn on curly vectors
resu@vcMinDistanceF = 0.03 ; thin the vectors
; resu@FieldTitle = "Wind " ; overwrite Field Title
; resu@UnitLabel = "m/s, vector"
; resu@vcRefAnnoOn = False ; show ref vector size value
resu@vcRefAnnoOn = True ; show ref vector size value
resu@vcMinAnnoOn = False
; resu@gsnLeftString = "Wind (vector)"
; plot vector every 2 grid point
jump = 1
vector =
wrf_vector(f0,wks,uzoom1(it,::jump,::jump),vzoom1(it,::jump,::jump),resu)
; Plot OLR in grey scale for the full domain
; plot = wrf_map_overlays(f0,wks,(/contour/),pltres,mpres)
; Now zoom into the box area and plot just for this area
plot = wrf_map_overlays(f0,wks,(/vector,contour1/),pltres,mpres)
; plot =
wrf_map_overlays(f0,wks,(/vector,contour1,contour2/),pltres,mpres)
; plot a box on the output
lats = (/ 46.0, 50.0 /)
lons = (/ -87.0, -83.0 /)
loc = wrf_user_ll_to_ij(f0, lons, lats, True)
lnres = True
lnres@gsLineColor = "NavyBlue"
lnres@gsLineThicknessF = 3.0
gsn_polyline(wks,plot,(/lons(0),lons(1)/),(/lats(0),lats(0)/),lnres)
gsn_polyline(wks,plot,(/lons(0),lons(1)/),(/lats(1),lats(1)/),lnres)
gsn_polyline(wks,plot,(/lons(0),lons(0)/),(/lats(0),lats(1)/),lnres)
gsn_polyline(wks,plot,(/lons(1),lons(1)/),(/lats(0),lats(1)/),lnres)
; plot a marker at the obs site
obsloc = wrf_user_ij_to_ll(f0, Ix, Iy, True) ; find lon and lat
print("lon/lat location is: " + obsloc)
lonobs = obsloc(0)
latobs = obsloc(1)
print("lon= " + lonobs + " lat= "+latobs)
gsres = True
gsres@tfPolyDrawOrder = "PreDraw"
; gsres@tfPolyDrawOrder = "PostDraw"
gsres@gsMarkerIndex = 16
gsres@gsMarkerColor = (/"black"/)
gsres@gsMarkerThicknessF = 10.7
; gsn_add_polymarker(wks,plot,lonobs,latobs,gsres)
gsn_polymarker(wks,plot,lonobs,latobs,gsres)
; put the plots together
resP = True
resP@gsnDraw = True ; draw
resP@gsnFrame = True ; advance frame
resP@gsnMaximize = True ; maximize image
resP@gsnPanelLabelBar = True ; add common colorbar
resP@gsnPanelTop = 1.0
resP@gsnPanelBottom = 0.1 ; space for label bar
resP@pmLabelBarWidthF = 0.95 ; make thinner
resP@pmLabelBarHeightF = 0.06
resP@txString = "The Gain of T at "+dateh
resP@lbLabelAutoStride = True
resP@gsnPaperOrientation = "portrait"
resP@pmLabelBarOrthogonalPosF = 0.2 ; move labelbar position
resP@gsnPanelXWhiteSpacePercent = 3
resP@gsnPanelYWhiteSpacePercent = 0
resP@gsnPanelYF = (/-2,-2,0.54,0.54/)
gsn_panel(wks,plot,(/1,1/), resP)
end do
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Feb 9 10:19:10 2011
This archive was generated by hypermail 2.1.8 : Fri Feb 11 2011 - 16:11:42 MST