; Example series of plotting meteograms with WRF ARW model data ; First let's just get and plot t2 at a point ; Add some inrfto to the plot ; Add slp to the plot ; Add a time-Z plot above slp and t2 ; Use wrf_user_ll_to_ij to get the point of interest ; Clean up the plot ;*********************************************** load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;*********************************************** procedure change_xaxis(wks,plot,taus,times) local res2, vpx, vpy, vph, vpw, blank begin ; ; The times array has values like "15_00","15_03",..., "15_21","16_00". ; ; For every "xx_00" label, we want a major tickmark, and for ; the others we want minor tickmark. ; ;---First, add the major tick marks. ii = str_match_ind(times,"_00") values = taus(ii) labels = str_get_cols(times(ii),0,1) setvalues plot "tmXBMode" : "Explicit" "tmXBValues" : values "tmXBLabels" : labels "tmXBLabelFontHeightF" : 0.012 end setvalues ;--Remove these so we can use them again delete(ii) delete(values) delete(labels) ; ; Second, get the locations of where we want minor tickmarks. ; We have to add these by creating a blank plot with a new ; set of tickarks. ; First, retrieve the viewport coordinates. This will be ; used for the blank plot. ; getvalues plot "vpXF" : vpx "vpYF" : vpy "vpHeightF" : vph "vpWidthF" : vpw end getvalues stimes = str_get_cols(times,3,4) ii = ind(stimes.ne."00") values = taus(ii) labels = stimes(ii) res2 = True res2@vpXF = vpx res2@vpYF = vpy res2@vpHeightF = vph res2@vpWidthF = vpw res2@tmXBMode = "Explicit" res2@tmXBValues = values res2@tmXBLabels = labels res2@tmXBLabelFontHeightF = 0.01 ; Make these labels smaller. res2@tmXBMajorOutwardLengthF = 0.005 ; Make the lengths smaller res2@tmXBMajorLengthF = 0.005 ; res2@tmXBLabelDeltaF = 0.6 ; Move label away from tickmarks. res2@tmXBLabelFontColor = "Brown" res2@tmYROn = False ; Turn off right tickmarks. res2@tmXTOn = False ; Turn off top tickmarks. res2@tmYLOn = False ; Turn off left tickmarks. blank = gsn_blank_plot(wks,res2) ; Create a blank plot. overlay(plot,blank) ; Overlay on existing plot. return end begin ;*********************************************** a = addfile("wrfout_d01.nc","r") ;----------------------------------------------------------------------- ; Find the ij location for the point if interest lat = -22.721949 lon = -65.695454 llres = True llres@ReturnInt = True ; Return integer values locij = wrf_user_ll_to_ij(a, lon, lat, llres) locij = locij - 1 ; array pointers in NCL space locX = locij(0) locY = locij(1) print(locY) print(locX) lat2d = a->XLAT(0,:,:) printMinMax(lat2d,0) lon2d = a->XLONG(0,:,:) printMinMax(lon2d,0) ; get time information and strip out the day and hour times_in_file = a->Times dims = dimsizes(times_in_file) times = new(dims(0),string) taus = ispan(1,dims(0),1) ; integer do i=0,dims(0)-1 times(i) = chartostring(times_in_file(i,8:12)) end do wks = gsn_open_wks("png","meteo_abrapampa") ; open a workstation ;----------------------------------------------------------------------- gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") rain_res = True t2 = wrf_user_getvar(a,"T2",-1) ; get t2 for all times t2=t2-273.15 t2@description="Temp.(C)" td2 = wrf_user_getvar(a,"td2",-1) ; get td2 for all times td2@description="P.Rocio(C)" slp = wrf_user_getvar(a,"slp",-1) ; get slp for all times hr = wrf_user_getvar(a,"rh2",-1) ; relative humidity hr@description="Hum. (%)" td = wrf_user_getvar(a,"td",-1) ; get tc for all times uvmet = wrf_user_getvar(a,"uvmet",-1) ; get rotated u and v comp of wind tc = wrf_user_getvar(a,"tc",-1) ; get tc for all times rain_exp = wrf_user_getvar(a,"RAINNC",-1) rain_con = wrf_user_getvar(a,"RAINC",-1) rain = rain_exp + rain_con u10=wrf_user_getvar(a,"U10",-1) v10=wrf_user_getvar(a,"V10",-1) velocity= (/ sqrt(u10^2+v10^2) /) velocity=velocity * 3.6 ;print (rain) temp = dimsizes(rain) print (temp) nt = temp(0) ;print (nt) new_rain = rain(1:nt-1,:,:) - rain(0:nt-2,:,:) ;print (new_rain) ;printVarSummary(t2) ;printVarSummary(slp) ;printVarSummary(tc) ;printVarSummary(uvmet) printVarSummary(t2) printVarSummary(hr) ;printVarSummary(new_rain) umet = uvmet(0,:,:,:,:) vmet = uvmet(1,:,:,:,:) ;----------------------------------------------------------------------- arain_point = rain(:,locY,locX) rain_point = new_rain(:,locY,locX) t2_point = t2(:,locY,locX) ; extract a time series at a point td2_point = td2(:,locY,locX) slp_point = slp(:,locY,locX) velocity_point=velocity(:,locY,locX) hr_point = hr(:,locY,locX) tc_point = tc(:,:,locY,locX) u_point = umet(:,:,locY,locX) v_point = vmet(:,:,locY,locX) t_max=(max(t2_point))+1 t_min=(min(td2_point))-1 t_nz=floattoint((t_max-t_min)/2) slp_max=(max(slp_point))+1 slp_min=(min(slp_point))-1 slp_nz= floattoint((slp_max-slp_min)/2) ar_max=(max(arain_point))+1 r_max=(max(rain_point))+1 rain_nz=floattoint((r_max-0)/2) wnd_max=(max(velocity_point))+1 ; Swap the dimensions as we want to plot time on the X axis later tt = tc_point(bottom_top|:,Time|:) ugrid = u_point(bottom_top|:,Time|:) vgrid = v_point(bottom_top|:,Time|:) ;----------------------------------------------------------------------- res2D = True ; Set basic resources res2D@gsnDraw = False ; Don't draw individual plot. res2D@gsnFrame = False ; Don't advance frame. res2D@vpXF = 0.15 ; x location res2D@vpYF = 0.95 ; y location res2D@vpWidthF = 0.65 ; width res2D@vpHeightF = 0.28 ; height res2D@tiXAxisString = "Dia_Hora" res2D@tiXAxisFontHeightF = 0.015 res2D@trYMaxF = 3 res2D@tmXBLabelJust = "CenterCenter" res2D@tmXBLabelFontHeightF = .015 res2D@tmLabelAutoStride = True tt_res = res2D tt_res@sfXArray = taus ;tt_res@gsnSpreadColors = True ; use full range of colors tt_res@cnFillOn = True ; turns on color fill tt_res@cnLevels = (/ -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45/) tt_res@cnFillColors = (/ 10, 16, 31, 40, 0, 55, 64, 80, 95, 110, 135, 159, 169, 191, 177/) tt_res@cnLevelSelectionMode = "ExplicitLevels" ; set levels manually tt_res@cnMinLevelValF = -20. tt_res@cnMaxLevelValF = 45. tt_res@cnLevelSpacingF = 5 tt_res@cnLinesOn = False tt_res@cnLineLabelsOn = False tt_res@cnInfoLabelOn = False tt_res@pmLabelBarDisplayMode = "Always" ; Add a label bar tt_res@pmLabelBarSide = "right" tt_res@pmLabelBarOrthogonalPosF = -0.05 tt_res@pmLabelBarParallelPosF = 0.04 tt_res@lbAutoManage = False tt_res@lbLabelAutoStride = True tt_res@lbOrientation = "vertical" tt_res@lbPerimOn = False tt_res@lbJustification = "BottomLeft" tt_res@lbLeftMarginF = 0.09 tt_res@lbTopMarginF = 0.09 tt_res@lbLabelFontHeightF = 0.012 tt_res@lbBoxLinesOn = False tt_res@tiYAxisString = "Altura (km)" tt_res@tiYAxisFontHeightF =0.015 tt_res@lbTitleOn =True tt_res@lbTitleString ="Temperatura" tt_res@lbTitleAngleF =0 tt_res@lbTitleJust ="CenterCenter" tt_res@lbTitlePosition = "right" tt_res@tiMainFontHeightF = 0.012 tt_res@tiMainString = "JUJUY: Perico" uv_res = res2D uv_res@vfXArray = taus uv_res@vcRefAnnoOn = False ; turns off the ref vector uv_res@vcRefLengthF = 0.090 ; set length of ref vector uv_res@vcGlyphStyle = "LineArrow" ; turn on wind barbs si quiero barbas cambio por WindBarb ; LineArrow flecha ;----------------------------------------------------------------------- res1D = True ; Set basic resources both will use res1D@vpXF = 0.15 ; The left side of the box location res1D@vpWidthF = 0.75 ; The Width of the plot box res1D@vpHeightF = 0.10 ; The height of the plot box res1D@xyLineThicknesses = 2 ; increase line thickness res1D@gsnDraw = False ; Don't draw individual plot. res1D@gsnFrame = False ; Don't advance frame. res1D@tmYLLabelFontHeightF =0.01 res1D@trXMinF = min(taus) res1D@trXMaxF = max(taus) res1D@tiYAxisFontHeightF=0.015 res1Dd = True res1Dd@gsLineDashPattern = 4 res1Dd@xyDashPattern = 4 res1Dd@xyDashPatterns = (/2,5,12/) res1Dd@trXMinF = min(taus) res1Dd@trXMaxF = max(taus) res1Dd@tiYAxisFontHeightF=0.015 slp_res = res1D slp_res@vpYF = 0.59 ; The top side of the plot box loc slp_res@xyLineColor = "red" ; set line color slp_res@tiYAxisString = "Pnm (hPa)" ; set y-axis string slp_res@tiYAxisFontColor = "red" hr_res = res1Dd hr_res@vpYF = 0.59 ; The top side of the plot box loc hr_res@xyLineColor = "blue" ; set line color hr_res@trYMinF = 0 ; min value on y-axis hr_res@trYMaxF = 100 ; min value on y-axis hr_res@tiYAxisFontColor = "blue" t2_res = res1D t2_res@vpYF = 0.95 ; The top side of the plot box loc t2_res@vpHeightF = 0.28 ; height t2_res@xyLineColor = "blue" ; set line color t2_res@trYMinF = t_min ; min value on y-axis t2_res@trYMaxF = t_max ; min value on y-axis t2_res@tmYLLabels = sprintf("%.0f",fspan(t_min,t_max,t_nz)) t2_res@tiYAxisFontColor = "blue" t2_res@tiMainFontHeightF = 0.015 t2_res@tiMainString = "JUJUY: Abrapampa" td2_res = res1Dd td2_res@vpYF = 0.44 ; The top side of the plot box loc td2_res@xyLineColor = "green" ; set line color td2_res@trYMinF = t_min ; min value on y-axis td2_res@trYMaxF = t_max ; min value on y-axis td2_res@tiYAxisFontColor = "green" ;----------------------------------------------------------------------- rain_res = res1D rain_res@vpYF = 0.59 ; The top side of the plot box location rain_res@vpHeightF = 0.28 ; height rain_res@tiXAxisString = "" ; X axis label. rain_res@tiYAxisString = "Prec. (mm)" ; Y axis label. rain_res@trYMinF = 0.0 ; min value on y-axis rain_res@trYMaxF = r_max ; min value on y-axis rain_res@xyLineThicknesses = 2 rain_res@gsnDraw = False ; Don't draw individual plot. rain_res@gsnFrame = False ; Don't advance frame. rain_res@gsnYRefLine = 0.0 ; create a reference line rain_res@gsnAboveYRefLineColor = "green" ; above ref line fill green rain_res@gsnXYBarChart = True ; turn on bar chart rain_res@tiYAxisFontColor = "green" arain_res = res1Dd arain_res@vpYF = 0.59 ; The top side of the plot box loc arain_res@vpHeightF = 0.28 ; height arain_res@xyLineColor = "blue" ; set line color arain_res@trYMinF = 0 ; min value on y-axis arain_res@trYMaxF = ar_max ; min value on y-axis arain_res@tiYAxisString = "P.Acum (mm)" ; Y axis label. arain_res@tiYAxisFontColor = "blue" wnd_res = res1D wnd_res@vpYF = 0.14 ; The top side of the plot box loc wnd_res@xyLineColor = "gold" ; set line color wnd_res@trYMinF = 0 ; min value on y-axis wnd_res@trYMaxF = wnd_max ; min value on y-axis wnd_res@tiYAxisString = "Viento 10m(km/h)" ; set y-axis string wnd_res@tiYAxisFontColor = "gold" ;----------------------------------------------------------------------- ; ttfill = gsn_contour(wks,tt,tt_res) ; change_xaxis(wks,ttfill,taus,times) ; windlayer = gsn_vector(wks,ugrid,vgrid,uv_res) ; overlay(ttfill,windlayer) ;wnd_plot = gsn_csm_xy(wks, taus, velocity_point, wnd_res) ;slp_plot = gsn_csm_xy2(wks,taus,slp_point,hr_point,slp_res,hr_res) t2_plot = gsn_csm_xy(wks,taus,t2_point,t2_res) rain_plot = gsn_csm_xy2(wks,taus,rain_point, arain_point,rain_res,arain_res) ; change_xaxis(wks,wnd_plot,taus,times) ; change_xaxis(wks,slp_plot,taus,times) change_xaxis(wks,t2_plot,taus,times) change_xaxis(wks,rain_plot,taus,times) ; draw(ttfill) ; draw(slp_plot) draw(t2_plot) draw(rain_plot) ;draw(wnd_plot) frame(wks) ;----------------------------------------------------------------------- end