;****************************************************** ; plotetaruc.ncl ;******************************************************* 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/csm/contributed.ncl" begin ;************************************************ ; open file and read in data ;************************************************ diri = "~/grib/" fili = "ruc2.t12z.pgrb20f06.grb" f = addfile (diri+fili, "r") names = getfilevarnames(f) ; print(names) ; printVarSummary(names) atts = getfilevaratts(f,names(0)) ; dims = getfilevardims(f,names(0)) ; print(atts) ; print(dims) ; get vars w = f->V_VEL_252_ISBL u = f->U_GRD_252_ISBL v = f->V_GRD_252_ISBL z = f->HGT_252_ISBL zeta = f->ABS_V_252_ISBL t = f->TMP_252_ISBL t_sfc = f->TMP_252_SFC u_sfc = f->U_GRD_252_HTGL v_sfc = f->V_GRD_252_HTGL pmsl = f->MSLMA_252_MSL rh = f->R_H_252_ISBL spd = (u^2.+v^2)^0.5 prate= f->PRATE_252_SFC pres = f->lv_ISBL2 lat2d = f->gridlat_252 lon2d = f->gridlon_252 time=w@initial_time print(time) ;printVarSummary(lat) dims=dimsizes(w) ; print(dims) nlon=dims(2) nlat=dims(1) ; print(nlon) ; print(nlat) ; print(var) ; print(pres(12)) ; printVarSummary(var) ; printVarSummary(lon2d) ; printVarSummary(lat2d) ;************************* ;common plotting resources ;************************* res = True ; plot mods desired res@gsnDraw = False res@gsnFrame = False res@cnLinesOn = False res@cnFillOn = True ; color plot desired res@cnLineLabelsOn = False ; turn off contour lines res@gsnAddCyclic = False ; !!!!! any plot of data that is on a native grid, must use the "corners" ; method of zooming in on map. ;**** ZOOMED IN res@mpLimitMode = "Corners" ; choose range of map res@mpLeftCornerLatF = 35. res@mpLeftCornerLonF = -105. res@mpRightCornerLatF = 47. res@mpRightCornerLonF = -80. ;**** ZOOMED OUT res@mpLimitMode = "Corners" ; choose range of map res@mpLeftCornerLatF = lat2d(0,0) res@mpLeftCornerLonF = lon2d(0,0) res@mpRightCornerLatF = lat2d(nlat-1,nlon-1) res@mpRightCornerLonF = lon2d(nlat-1,nlon-1) ; ; The following 4 pieces of information are REQUIRED to properly display ; data on a native lambert conformal grid. This data should be specified ; somewhere in the model itself. ; WARNING: our local RCM users could not provide us with this information, ; so this is our best guess as to the correct numbers. Use at your own risk. res@mpProjection = "LambertConformal" res@mpLambertParallel1F = lon2d@Latin1 res@mpLambertParallel2F = lon2d@Latin2 res@mpLambertMeridianF = lon2d@Lov ; usually, when data is placed onto a map, it is TRANSFORMED to the specified ; projection. Since this model is already on a native lambert conformal grid, ; we want to turn OFF the tranformation. res@tfDoNDCOverlay = True res@mpGeophysicalLineColor = "red" ; color of continental outlines res@mpPerimOn = True ; draw box around map res@mpGridLineDashPattern = 2 ; lat/lon lines as dashed res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; add state boundaries res@mpUSStateLineColor = "red" ; make them red res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks res@cnLineLabelInterval = 1 res@cnLabelMasking = True res@gsnSpreadColors=True res@gsnSpreadColorStart = 3 res@gsnSpreadColorEnd = 19 wks = gsn_open_wks("ps","ruc_300_zoomed") ; open a workstation ;plot=new(4,graphic) gsn_define_colormap(wks,"gui_default") ; choose colormap ;************************************************ ; create plot: upper left ;************************************************ res1=res res1@gsnMaximize = True res1@tiMainString=time res1@gsnLeftString="300hPa height(dm), wind speed(m/s) " res1@cnLevelSelectionMode = "ManualLevels" ; manual levels res1@cnMinLevelValF = 10. ; min level res1@cnMaxLevelValF = 100. ; max level res1@cnLevelSpacingF = 10 ; contour spacing res2=True res2@tfDoNDCOverlay = True res2@cnInfoLabelOn = False res2@cnMinLevelValF = 900.+12*12 ; set min contour level res2@cnMaxLevelValF = 900.-12*12 ; set max contour level res2@cnLevelSpacingF = 12. ; set contour spacing res2@cnLineLabelsOn = True res2@cnLineLabelInterval = 1 res2@cnLabelMasking = True ; res2@cnHighLabelsOn = True; label highs ; res2@cnHighLabelString = "H"; highs' label ; res2@cnHighLabelFontHeightF = 0.03; larger H labels ; res2@cnHighLabelFont = "helvetica-bold"; H labels font ; res2@cnHighLabelBackgroundColor = "Transparent"; no box ; res2@cnLowLabelsOn = True; label lows ; res2@cnLowLabelString = "L"; lows' label ; res2@cnLowLabelFontHeightF = 0.03; larger L labels ; res2@cnLowLabelFont = "helvetica-bold"; L labels font ; res2@cnLowLabelBackgroundColor = "Transparent"; no box ; res2@cnConpackParams = (/ "HLX:6, HLY:6" /) index300=ind(pres.eq.300.) plot1 = gsn_csm_contour_map_overlay(wks,spd(index300,:,:),z(index300,:,:)/10.,res1,res2) ; Draw contours over a map. ; draw(plot1) ; frame(wks) getvalues plot1 "mpLeftNPCF" : lnpc "mpBottomNPCF" : bnpc "mpTopNPCF" : tnpc "mpRightNPCF" : rnpc end getvalues print("NPC (before zooming) bottom: " + bnpc + " left: " + lnpc + " top: " + tnpc + " right: " + rnpc) res1@mpLeftCornerLatF = 35. res1@mpLeftCornerLonF = -105. res1@mpRightCornerLatF = 47. res1@mpRightCornerLonF = -80. plot1 = gsn_csm_contour_map_overlay(wks,spd(index300,:,:),z(index300,:,:)/10.,res1,res2) ; Draw contours over a map. getvalues plot1 "mpLeftNPCF" : lnpc1 "mpBottomNPCF" : bnpc1 "mpTopNPCF" : tnpc1 "mpRightNPCF" : rnpc1 end getvalues print("NPC (after zooming) bottom: " + bnpc1 + " left: " + lnpc1 + " top: " + tnpc1 + " right: " + rnpc1) width = rnpc - lnpc lf = (lnpc1 - lnpc) / width rf = (rnpc1 - lnpc) / width height = tnpc - bnpc bf = (bnpc1 - bnpc) / height tf = (tnpc1 - bnpc) / height print("Fractional distance from bottom left, left: " + lf + " right: " + rf + " bottom: " + bf + " top: " + tf) dims = dimsizes(spd) xgrid = dims(2) ygrid = dims(1) xl = floattoint(lf * xgrid) xr = floattoint(rf * xgrid) yl = floattoint(bf * ygrid) yr = floattoint(tf * ygrid) print("Grid subscripting indexes, x left: " + xl + " x right: " + xr + " y left: " + yl + " yright: " + yr) print("Coordinates at these subscripts, bottom left lat/lon: " \ + f->gridlat_252(yl,xl) + "," + f->gridlon_252(yl,xl) + \ " top right lat/lon: " + f->gridlat_252(yr,xr) + "," + f->gridlon_252(yr,xr)) plot1 = gsn_csm_contour_map_overlay(wks,spd(index300,yl:yr,xl:xr),z(index300,yl:yr,xl:xr)/10.,res1,res2) ; Draw contours over a map. draw(plot1) frame(wks) end