; 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 = "maxwind1day.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 (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,"BlAqGrYeOrReVi200") maxvalue = max(sitedata) minvalue = min(floor(sitedata)) minvalue = minvalue - mod(minvalue,5) maxvalue = ceil((5-mod(maxvalue,5))+maxvalue) diff = toint((maxvalue-minvalue))+2 if diff.gt.99 colors = round(fspan(2,199,diff),3) res@cnLevelSpacingF = 1.0 ;-- contour level spacing res@lbLabelStride = 10.0 end if if (diff.lt.99).and.(diff.ge.49) colors = round(fspan(2,199,2*diff),3) res@cnLevelSpacingF = 0.5 ;-- contour level spacing res@lbLabelStride = 10.0 end if if diff.lt.49 colors = round(fspan(2,199,4*diff),3) res@cnLevelSpacingF = 0.25 ;-- contour level spacing res@lbLabelStride = 20.0 end if res@cnMinLevelValF = minvalue ;-- minimum contour level res@cnMaxLevelValF = 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 diff.ge.25 cres@cnLevels = ispan(-70,120,10) end if if diff.lt.25 cres@cnLevels = ispan(-70,120,5) 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 = "Maximum Wind Gust (mph)" gsn_text_ndc(wks,label,.32,.35 ,txres) txres@txFontHeightF = .022 label = "On " + systemfunc("date -u -d -1days +%A") + ", " + systemfunc("date -u -d -1days +%F") gsn_text_ndc(wks,label,.32,.29 ,txres) frame(wks) end