; load libraries, twice for good measure 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" ;load "/home/torn/ncl/wrf_map_info.ncl" ;load "/home/torn/ncl/array_stat.ncl" ;load "/home/torn/ncl/image_conv.ncl" load "shapefile_utils.ncl" begin ;---Read in site info sitefile = "precip2day.dat" sitedatafile = asciiread(sitefile,-1,"string") delim = "," sitelats = tofloat(str_get_field(sitedatafile,1,delim)) sitelons = tofloat(str_get_field(sitedatafile,2,delim)) sitedata = tofloat(str_get_field(sitedatafile,3,delim)) ;;;; this section gets mask for risk/colorbar purposes in only the domain sname = "USA_adm1.shp" maskres = True maskres@return_mask = True maskres@shape_var = "NAME_1" maskres@shape_names = "New York" maskres@debug = True f=addfile("MRMS_GaugeCorr_QPE_48H_latest.grib2","r") newdata = f->VAR_209_6_35_P0_L102_GLL0({40.3:45.3},{280:288.3}) newdata&lon_0=newdata&lon_0-360 maxcheck = newdata(::5,::5) nymask := shapefile_mask_data(maxcheck,sname,maskres) maxcheck = maxcheck*nymask*0.0393701 newdata =newdata*0.0393701 delete_VarAtts(newdata,"units") delete_VarAtts(newdata,"long_name") wks_type = "png" wks_type@wkWidth = 4000 wks_type@wkHeight = 3200 wks = gsn_open_wks(wks_type,"plot") res = True res@gsnAddCyclic = False res@mpProjection = "Mercator" res@mpProjection = "CylindricalEquidistant" res@mpLimitMode = "Corners" ;res@mpOutlineBoundarySets = "GeophysicalAndUSStates" res@mpLeftCornerLatF = 40.45 res@mpLeftCornerLonF = -80 res@mpRightCornerLatF = 45.2 res@mpRightCornerLonF = -71.8 res@cnFillOn = True res@cnFillMode = "RasterFill" res@cnRasterSmoothingOn = True res@cnLinesOn = False ;no contours on top of fill res@gsnMaximize = True res@gsnFrame = False res@tmXBBorderOn = True res@tmBorderLineColor = "white" res@tmXBMinorPerMajor = 0 res@tmXBOn = False res@tmXTBorderOn = True res@tmXBMajorLineColor= "white" res@tmXTOn = False res@tmYLBorderOn = True res@tmYLMajorLineColor= "white" res@tmYLMinorPerMajor = 0 res@tmYLOn = False res@tmYRBorderOn = True res@tmYRMajorLineColor= "white" res@tmYROn = False res@mpDataBaseVersion = "Ncarg4_1" res@mpDataSetName = "Earth..4" ; For counties res@mpFillAreaSpecifiers = (/"Land", "Land:North America:United States:Conterminous US:New York","Water"/) res@mpSpecifiedFillColors = (/"white","transparent", "white"/) res@mpInlandWaterFillColor = 238 res@mpLandFillColor = -1 ; transparent res@mpOutlineOn = True res@mpMaskOutlineSpecifiers = (/"Land"/) res@mpOutlineSpecifiers = (/ "Land:North America:United States:Conterminous US:New York:counties"/) res@mpDataBaseVersion = "MediumRes" ;-- better map resolution res@mpGeophysicalLineThicknessF = 18 res@mpUSStateLineThicknessF = 18 res@mpNationalLineThicknessF = 18 res@mpCountyLineThicknessF = 5 res@cnFillDrawOrder = "PreDraw" res@gsnDraw = False res@cnLevelSelectionMode = "ManualLevels" ;-- set contour levels manually res@mpOutlineDrawOrder = "Draw" gsn_define_colormap(wks,"wh-bl-gr-ye-re") maxvalue = where(max(maxcheck).gt.(.10),max(maxcheck)+.1*max(maxcheck),.10) minvalue = 0 ;minvalue = min(floor(sitedata)) ;minvalue = minvalue - mod(minvalue,5) ;maxvalue = ceil((5-mod(maxvalue,5))+maxvalue) if ceil(100*maxvalue).lt.10 colors := round(fspan(2,199,toint(maxvalue*100)+10),3) res@cnLevelSpacingF = 0.005 res@lbLabelStride = 1.0 end if if (ceil(100*maxvalue).lt.195).and.(ceil(100*maxvalue).ge.10) colors := round(fspan(2,199,toint(maxvalue*100)+5),3) res@cnLevelSpacingF = 0.01 res@lbLabelStride = 5.0 end if if (ceil(100*maxvalue).ge.195).and.(ceil(100*maxvalue).lt.500) colors := round(fspan(2,199,toint(maxvalue*50)+2),3) res@cnLevelSpacingF = 0.02 res@lbLabelStride = 10.0 end if if ceil(100*maxvalue).ge.395 colors := round(fspan(2,199,toint(maxvalue*40)+2),3) res@cnLevelSpacingF = 0.025 res@lbLabelStride = 10.0 end if if ceil(100*maxvalue).ge.595 colors := round(fspan(2,199,toint(maxvalue*30)+2),3) res@cnLevelSpacingF = 0.04 res@lbLabelStride = 25.0 end if if ceil(100*maxvalue).ge.795 colors := round(fspan(2,199,toint(maxvalue*20)+2),3) res@cnLevelSpacingF = 0.05 res@lbLabelStride = 20.0 end if colors := array_append_record(0,colors,0) res@cnMinLevelValF = 0 ;-- minimum contour level res@cnMaxLevelValF = .01*(ceil(100*maxvalue)) ;-- maximum contour level res@lbOrientation = "Vertical" res@cnFillColors = colors res@lbBoxLinesOn = False res@lbTitleString = " " res@lbTitleFontHeightF = 0.012 res@lbLabelFontHeightF = 0.01 res@pmLabelBarOrthogonalPosF = -0.20 res@pmLabelBarParallelPosF = 0.65 res@pmLabelBarHeightF = .3 plot = gsn_csm_contour_map(wks,newdata,res) ;;;;;;;;;;;; this part tells the state marker what color to be based on the data range := fspan(minvalue,maxvalue,toint((maxvalue-minvalue)/res@cnLevelSpacingF)) markercolors := new(dimsizes(sitedata),"integer") markercolors := new(dimsizes(sitedata),"integer") do n=0,dimsizes(sitedata)-1 tempvalues := where((sitedata(n)-range).lt.0,999,(sitedata(n)-range)) markercolors(n) = colors(minind(tempvalues)) end do markercolors=where(sitedata.gt.maxvalue,0,markercolors) markercolors=where(sitedata.lt.minvalue,1,markercolors) stationmarker := new(dimsizes(sitedata),graphic) blackmarkers := stationmarker mres = True ; marker resources for best track mres@gsMarkerIndex = 16 ; marker style (filled circle) mres@gsMarkerSizeF = 16.0 ; marker size mres@gsMarkerColor = "gray50" ; marker color mres@tfPolyDrawOrder= "PostDraw" blackmarkers := gsn_add_polymarker(wks,plot,sitelons,sitelats,mres) mres@gsMarkerSizeF = 25 keymarkers = gsn_add_polymarker(wks,plot,(/-79.8,-77/),(/40.65,40.65/),mres) markercolors=where(sitedata.gt.maxvalue,1,markercolors) do n=0,dimsizes(sitedata)-1,1 mres@gsMarkerSizeF = 10 ; marker size mres@gsMarkerColor := markercolors(n) ; marker color stationmarker(n) = gsn_add_polymarker(wks,plot,sitelons(n),sitelats(n),mres) end do mres@gsMarkerSizeF = 15 ; marker size mres@gsMarkerColor = round(avg(colors),3) plusmarkers1 = gsn_add_polymarker(wks,plot,-79.8,40.65,mres) mres@gsMarkerColor = 0 ; plusmarkers0 = gsn_add_polymarker(wks,plot,-79.8,41.2,mres) mres@gsMarkerColor = 1 plusmarkers2 = gsn_add_polymarker(wks,plot,-77,40.65,mres) txres = True txres@txFont = "helvetica-bold" txres@txFontHeightF = .013 ; Set the font height txres@txJust = "CenterLeft" plustext1 = gsn_add_text(wks,plot,"Equivalent Mesonet Observation",-79.65,40.65,txres) ;; add units text plustext0 = gsn_add_text(wks,plot,"Observation Above NY MRMS Max",-76.85,40.65,txres) ;; add units text ; plustext2 = gsn_add_text(wks,plot,"Observation Below Statewide Minimum",-79.6,40.8,txres) ;; add units text if any(markercolors.eq.1) txres@txFontHeightF = .0095 txres@txFontColor = "magenta" quote = inttochar(34) abovelats = sitelats(ind(markercolors.eq.1)) abovelons = sitelons(ind(markercolors.eq.1)) abovedata = sitedata(ind(markercolors.eq.1)) abovetext = gsn_add_text(wks,plot,.01*round((100*abovedata),3)+""+quote,abovelons+.07,abovelats-.07,txres) end if draw(plot) txres = True txres@txFontHeightF = .03 ; Set the font height txres@txJust = "CenterCenter" txres@txFontColor = "black" label = "Total Precipitation (in)" gsn_text_ndc(wks,label,.32,.36 ,txres) txres@txFontHeightF = .022 label = "Days Of: " + systemfunc("date -u -d -2days +%F") + " And " + systemfunc("date -u -d -1days +%F") gsn_text_ndc(wks,label,.32,.315 ,txres) label = "Fill Is MRMS Estimated Precipitation" gsn_text_ndc(wks,label,.32,.28 ,txres) frame(wks) end