; Plot the observed precipitation for the Feb 21 storm of the third OSE case in 2016 ; Uses radar and station data to plot the 24 hour and accumulated precipitation plots ; Using shape files ; Andrew Kren ; May 27, 2016 ;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; shape file example extensions ; nws_precip_1day_observed_20160221.dbf ; nws_precip_1day_observed_20160221.shx ; nws_precip_1day_observed_20160221.shp ; nws_precip_1day_observed_20160221.prj ;;;; begin obs_dir = "/scratch4/BMC/shout/Andrew.Kren/ENRR_obs_precip/" fils_dbf = systemfunc("csh -c 'cd " + obs_dir + " ; ls *dbf'") fils_shx = systemfunc("csh -c 'cd " + obs_dir + " ; ls *shx'") fils_shp = systemfunc("csh -c 'cd " + obs_dir + " ; ls *shp'") fils_prj = systemfunc("csh -c 'cd " + obs_dir + " ; ls *prj'") lat_low = 32 ; coordinates used to plot the verification region for CA lat_high = 46 lon_low = 231 lon_high = 244 lat_low_2 = 55 ; coordinates used to plot the verification region for AK lat_high_2 = 65 lon_low_2 = 193 lon_high_2 = 220 ; --------------------------- end of user's settings ----------------------------------------- do i=0,dimsizes(fils_dbf)-1 print("Reading Shapefile " + fils_dbf(i)) w = addfile(obs_dir+fils_shp(i),"r") id = w->ID hrapx = w->HRAPX hrapy = w->HRAPY lat = w->LAT lon = w->LON pcp = w->GLOBVALUE lon = lon + 360.0 printVarSummary(pcp) exit pcp!0="time" mslp_nr!1="lat_3" mslp_nr!2="lon_3" mslp_nr&lat_3 = 1.0*ispan(90,-90,1) mslp_nr&lon_3 = 1.0*ispan(0,359,1) mslp_nr&lat_3@units = "degrees_north" mslp_nr&lon_3@units = "degrees_east" mslp_nr@_FillValue=1.e20 delete(w) w = addfile(obs_dir+fils_prj(i),"r") lat = w->LATITUDE lon = w->LONGITUDE exit end do ;hourly_pcp = new((/dimsizes(fils),nlat,nlon/),float) ; 24 hr accumulated precip for the current day ; hourly_pcp!0="hour" ; hourly_pcp!1="lat" ; hourly_pcp!2="lon" ; hourly_pcp@_FillValue = -999 ;accum_pcp = new((/dimsizes(fils),nlat,nlon/),float) ; accumulated precipitation over each day ; accum_pcp!0="hour" ; accum_pcp!1="lat" ; accum_pcp!2="lon" ; accum_pcp@_FillValue = -999 if (region.eq."conus") then year = tointeger(str_get_cols(fils,17,20)) ; year, month, day, hour, minute, sec in UTC, same for hamsr month = tointeger(str_get_cols(fils,21,22)) day = tointeger(str_get_cols(fils,23,24)) else year = tointeger(str_get_cols(fils,14,17)) ; year, month, day, hour, minute, sec in UTC, same for hamsr month = tointeger(str_get_cols(fils,18,19)) day = tointeger(str_get_cols(fils,20,21)) end if do i=0,dimsizes(fils)-1 ; loop over all daily files print("Reading YMD " + tostring(year(i))+tostring(month(i))+tostring(day(i))) print("Reading file: " + fils(i)) w = addfile(obs_dir+fils(i),"r") lat_corners = w->lat lon_corners = w->lon true_lat = w->true_lat true_lon = w->true_lon factor = 0.0003937008 pcp = factor*tofloat(w->amountofprecip) ; to inches hourly_pcp(i,:,:) = (/pcp/) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; beginning of plot section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; outfile = "pcp_" + region + "_" + tostring(year(i)) + tostring(month(i)) + tostring(day(i)) wks = gsn_open_wks("eps",outfile) gsn_define_colormap(wks,"WhBlGrYeRe") res = True res@mpProjection = "Stereographic" res@mpLimitMode = "Corners" res@mpLeftCornerLatF = lat_corners(0) res@mpLeftCornerLonF = lon_corners(0) res@mpRightCornerLatF = lat_corners(2) res@mpRightCornerLonF = lon_corners(2) res@mpEllipticalBoundary = True res@mpRelativeCenterLon = True ; set a center lon res@mpCenterLonF = true_lon ; center lon res@mpRelativeCenterLat = True ; set a center lat res@mpCenterLatF = true_lat ; center lat res@tfDoNDCOverlay = True res@cnLinesOn = False res@gsnMaximize = True res@cnFillOn = True ; color plot desired res@cnLineLabelsOn = False ; turn off contour lines res@tiMainString = "24-hr Accumulated Precipitation: " + tostring(month(i)) + "/" + tostring(day(i)) + "/" + tostring(year(i)) ;res@mpMinLatF = 15 ; specify min lat ;res@mpMaxLatF = 70 ;res@mpMinLonF = 150 ;res@mpMaxLonF = 270 ;res@mpCenterLonF = 220 res@gsnRightString = "" res@gsnLeftString = "" res@mpGridAndLimbOn = True res@mpGridLatSpacingF = 10.0 res@mpGridLonSpacingF = 10.0 res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ;res@mpOutlineBoundarySets = "National" ; turn on country boundaries res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries res@mpUSStateLineThicknessF = 2.5 res@gsnFrame = False ; don't advance frame res@gsnDraw = False ; don't draw plot res@cnFillOpacityF = 0.9 res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/ 0.01,0.05,0.10,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4/) res@lbLabelStrings = (/"0.01","0.05","0.10", "0.25", "0.5", "0.75", "1", "1.5", "2","2.5","3","3.5","4"/) res@cnExplicitLabelBarLabelsOn = True res@lbTitleOn = True res@lbTitleString = "Precipitation (in)" ; title string res@lbTitlePosition = "Top" ; title position res@lbLabelFontHeightF =.022 ; make labels larger res@lbTitleFontHeightF= .015 ; make title smaller res@pmLabelBarOrthogonalPosF = .10 ; move whole thing down plot = gsn_csm_contour_map(wks,hourly_pcp(i,:,:),res) gres = True gres@gsLineColor = "Red" gres@gsLineThicknessF = 5.0 lat2 = (/ lat_low, lat_low, lat_high, lat_high, lat_low/) lon2 = (/lon_low, lon_high, lon_high, lon_low, lon_low/) dum = gsn_add_polyline(wks,plot,lon2,lat2,gres) draw(plot) frame(wks) ; advance frame delete(res) delete([/lat_corners,lon_corners,true_lat,true_lon,pcp/]) end do print("Calculating accumulated precipitation") ntstrt = 0 ; counter variables for accumulated precipitation ntlast = 0 do t=0,dimsizes(fils)-1 accum_pcp(t,:,:) = dim_sum_n(hourly_pcp(ntstrt:ntlast,:,:),0) ntlast = ntlast + 1 ; increment for accumulated precipitation end do exit print("Making accumulated pcp plots...") ;;;;;;;;;;;;;;;;;;; ; Plot accumulated precipitation ;;;;;;;;;;;;;;;;;;; ;do i=0,fc_files-1 ;end do print("Program ended successfully") end