; 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 = "weeklylowanomaly.dat" sitedata = asciiread(sitefile,-1,"string") delim = "," sitelats = tofloat(str_get_field(sitedata,1,delim)) sitelons = tofloat(str_get_field(sitedata,2,delim)) siteprecip= tofloat(str_get_field(sitedata,3,delim)) 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,siteprecip,newlats,newlons) newdata2 = dsgrid2(sitelats,sitelons,siteprecip,newlats,newlons) newdata = natgrid(sitelats,sitelons,siteprecip,newlats,newlons) ;newdata = obj_anal_ic(sitelats,sitelons,siteprecip,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,"weeklylowanomalymap") 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"/) res@mpDataBaseVersion = "MediumRes" ;-- better map resolution res@mpGeophysicalLineThicknessF = 4 res@mpUSStateLineThicknessF = 4 res@mpNationalLineThicknessF = 4 res@cnFillDrawOrder = "PreDraw" res@gsnDraw = False res@cnLevelSelectionMode = "ManualLevels" ;-- set contour levels manually ;gsn_define_colormap(wks,"MPL_YlGn") colors = round(fspan(27,227,162),3) res@cnMinLevelValF =-10 ; minvalue ;-- minimum contour level res@cnMaxLevelValF = 10 ; maxvalue ;-- maximum contour level res@cnLevelSpacingF = 0.125 ;-- contour level spacing res@lbLabelStride = 16 res@lbOrientation = "Vertical" res@cnFillColors = colors res@lbBoxLinesOn = False res@lbTitleString = "Anomaly (F)" res@lbTitleFontHeightF = 0.015 res@lbLabelFontHeightF = 0.01 res@pmLabelBarOrthogonalPosF = -0.18 res@pmLabelBarParallelPosF = 0.65 res@pmLabelBarHeightF = .4 hghtfile = "topo15.dat" hghtdatafile = asciiread(hghtfile,-1,"string") hghtlats = tofloat(str_get_field(hghtdatafile,3,delim)) hghtlons = tofloat(str_get_field(hghtdatafile,4,delim)) hght = tofloat(str_get_field(hghtdatafile,5,delim))*3.28084 ;newdata = cssgrid(sitelats,sitelons,siteprecip,newlats,newlons) hghtdata2 = dsgrid2(hghtlats,hghtlons,hght,newlats,newlons) hghtdata = natgrid(hghtlats,hghtlons,hght,newlats,newlons) ;newdata = obj_anal_ic(sitelats,sitelons,siteprecip,newlats,newlons,(/2,1,.5/),False) hghtdata = (hghtdata+hghtdata2)/2 hghtdata!0="lat" hghtdata!1="lon" hghtdata&lon = newlons hghtdata&lat = newlats cres = True cres@cnLevelSelectionMode = "ExplicitLevels" cres@cnLevels = (/500,750,1000,1250,1500,1750,2000,2250,2500/) ;cres@cnLineThicknesses = (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21/) cres@cnLineThicknesses = (/5,6,7,8,9,10,11,12,13,14,15,20,22/) cres@cnInfoLabelOn = False cres@cnLineColors = (/"gray2","gray10","gray18","gray26","gray34","gray42","gray50","gray58","gray64","gray72"/) cres@cnLineLabelFontColors = "red" cres@cnMonoLineLabelFontColor = True cres@cnMonoLineThickness = False cres@cnMonoLineColor = False cres@cnMonoLevelFlag = False cres@cnLineLabelsOn = True cres@cnLineLabelBackgroundColor = -1 ; transparent cres@cnLineLabelPlacementMode = "Constant" cres@cnLineLabelInterval = 0 cres@gsnAddCyclic = False cres@cnLineDashSegLenF = .3 cres@cnLineDrawOrder = "PreDraw" plot = gsn_csm_contour_map_overlay(wks,newdata,hghtdata,res,cres) ;mres = True ; marker resources for best track ;mres@gsMarkerIndex = 16 ; marker style (filled circle) ;mres@gsMarkerSizeF = 5.0 ; marker size ;mres@gsMarkerColor = "gray" ; marker color ;markers = gsn_add_polymarker(wks,plot,sitelons,sitelats,mres) draw(plot) txres = True txres@txFontHeightF = .015 ; Set the font height txres@txJust = "CenterLeft" label = "Weekly Low Temperature Anomaly (F, fill) and Elevation" gsn_text_ndc(wks,label,.03,.34 ,txres) label = "Time Period: "+systemfunc("date -u -d -7days +%a")+"., "+systemfunc("date -u -d -7days +%F")+" through "+ systemfunc("date -u -d -1days +%a")+"., "+systemfunc("date -u -d -1days +%F") gsn_text_ndc(wks,label,.03,.30 ,txres) label = (dimsizes(sitelats))+" Stations Were Used To Make This Plot" gsn_text_ndc(wks,label,.03,.26 ,txres) frame(wks) end