; 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 "$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 "wrf_dart_util.ncl" begin ;---Read in site info sitefile = "snow2day.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 part adds fake stations around the periphery of NY so that the map is smoother near the border size = dimsizes(sitelats)-1 do i = 1,size if sitelats(i).gt.44.87 sitelats := array_append_record(sitelats,sitelats(i)+1,0) sitelons := array_append_record(sitelons,sitelons(i),0) sitedata := array_append_record(sitedata,sitedata(i),0) end if if (sitelats(i).lt.42.1) .and. (sitelons(i).lt.(-75.9)) sitelats := array_append_record(sitelats,sitelats(i)-1,0) sitelons := array_append_record(sitelons,sitelons(i),0) sitedata := array_append_record(sitedata,sitedata(i),0) end if if (sitelats(i).gt.41.5) .and. (sitelons(i).gt.(-73.6)) sitelats := array_append_record(sitelats,sitelats(i),0) sitelons := array_append_record(sitelons,sitelons(i)+1,0) sitedata := array_append_record(sitedata,sitedata(i),0) end if if (sitelons(i).lt.(-78.74)) sitelats := array_append_record(sitelats,sitelats(i),0) sitelons := array_append_record(sitelons,sitelons(i)-1,0) sitedata := array_append_record(sitedata,sitedata(i),0) end if end do lonmin = -81 ;round(min(sitelons),3)-2 lonmax = -71 ;round(max(sitelons),3)+2 latmin = 38 ;round(min(sitelats),3)-2 latmax = 48 ;round(max(sitelats),3)+2 newlons = fspan(lonmin,lonmax,400) newlats = fspan(latmin,latmax,400) latsize=(dimsizes(newlats)) lonsize=(dimsizes(newlons)) newlons@units="degrees_east" newlats@units="degrees_north" ;newdata = cssgrid(sitelats,sitelons,sitedata,newlats,newlons) ;newdata2 = dsgrid2(sitelats,sitelons,sitedata,newlats,newlons) newdata = natgrid(sitelats,sitelons,sitedata,newlats,newlons) ;newdata = obj_anal_ic(sitelats,sitelons,sitedata,newlats,newlons,(/2,1,.5/),False) ;newdata = (newdata+newdata2)/2 newdata!0="lat" newdata!1="lon" newdata&lon = newlons newdata&lat = newlats wks_type = "png" wks_type@wkWidth = 2000 wks_type@wkHeight = 1600 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 = 9 res@mpUSStateLineThicknessF = 9 res@mpNationalLineThicknessF = 9 res@mpCountyLineThicknessF = 2.5 res@cnFillDrawOrder = "PreDraw" res@gsnDraw = False res@cnLevelSelectionMode = "ManualLevels" ;-- set contour levels manually gsn_define_colormap(wks,"wh-bl-gr-ye-re") maxvalue = max(sitedata) minvalue = 0 ;minvalue = min(floor(sitedata)) ;minvalue = minvalue - mod(minvalue,5) ;maxvalue = ceil((5-mod(maxvalue,5))+maxvalue) if ceil(maxvalue).le.4 colors := round(fspan(2,199,toint(maxvalue*15)+5),3) res@cnMinLevelValF = .5 ;-- minimum contour level res@cnLevelSpacingF = 0.1 res@lbLabelStride = 5.0 end if if (ceil(maxvalue).lt.10).and.(ceil(maxvalue).ge.4) colors := round(fspan(2,199,toint(maxvalue*4)+10),3) res@cnMinLevelValF = 1 ;-- minimum contour level res@cnLevelSpacingF = 0.2 res@lbLabelStride = 5.0 end if if (ceil(maxvalue).ge.10).and.(ceil(maxvalue).lt.25) colors := round(fspan(2,199,toint(4*maxvalue)+2),3) res@cnMinLevelValF = 1 ;-- minimum contour level res@cnLevelSpacingF = 0.25 res@lbLabelStride = 16.0 end if if ceil(maxvalue).ge.25 colors := round(fspan(2,199,toint(2*maxvalue)+2),3) res@cnMinLevelValF = 2 ;-- minimum contour level res@cnLevelSpacingF = 0.5 res@lbLabelStride = 8.0 end if ;colors := array_append_record(0,colors,0) res@cnMaxLevelValF = (ceil(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 cres = True cres@cnLevelSelectionMode = "ExplicitLevels" if (ceil(maxvalue)).lt.4 cres@cnLevels = .5*ispan(1,10,1) end if if ((ceil(maxvalue)).ge.4) .and. ((ceil(maxvalue)).le.10) cres@cnLevels := ispan(1,11,1) end if if (ceil(maxvalue).ge.10).and.(ceil(maxvalue).lt.25) cres@cnLevels := ispan(3,27,3) end if if (ceil(maxvalue)).ge.25 .and.(ceil(maxvalue).lt.50) cres@cnLevels := ispan(5,800,5) end if if (ceil(maxvalue)).ge.50 cres@cnLevels := ispan(10,800,10) end if cres@cnLineThicknessF = 12.5 cres@cnInfoLabelOn = False cres@cnLineColors = "red" cres@cnLineLabelFontColors = "red" cres@cnMonoLineLabelFontColor = True cres@cnMonoLineThickness = True cres@cnMonoLineColor = True cres@cnMonoLevelFlag = False cres@cnLineLabelsOn = True cres@cnLineLabelBackgroundColor = -1 ; transparent cres@cnLineLabelPlacementMode = "Constant" cres@cnLineLabelInterval = 1 cres@gsnAddCyclic = False cres@cnLineDashSegLenF = .1 cres@cnLineDrawOrder = "PreDraw" plot = gsn_csm_contour_map_overlay(wks,newdata,newdata,res,cres) draw(plot) txres = True txres@txFontHeightF = .03 ; Set the font height txres@txJust = "CenterCenter" label = "Total Snowfall (in)" gsn_text_ndc(wks,label,.32,.35 ,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,.29 ,txres) frame(wks) end