; 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  = "precip1day.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_24H_latest.grib2","r")
newdata = f->GaugeCorrQPE24H_P0_L102_GLL0({40.3:45.3},{280:288.3})
newdata&lon_0=newdata&lon_0-360
maxcheck = newdata(::6,::6)
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 = "On " + systemfunc("date -u -d -1days +%A") + ", " + 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