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/csm/contributed.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "./WRFUserARW.ncl"
load "./shapefile_mask_data.ncl"

begin

; Open up the workstation so we can make a plot
wks  = gsn_open_wks ("X11","shape_mask")

; Open up the T2 max/min file that we just made
a = addfile("t2_maxmin_daily_jja_1990_rtty_1990.nc","r")
; Read in t2max, lat and lon
t2max = a->T2MAX(0,:,:)
t2max@lat2d = a->XLAT
t2max@lon2d = a->XLONG

; What is our shapefile name?
DIRs = "data/ShapeFiles_US/"
shp_filename = DIRs + "cb_2014_us_state_20m.shp"

; Let's find out what is in the shapefile!
f = addfile(shp_filename,"r")
print(f)
id = f->NAME
print(id)

; Setting the options to rad the shapefile 
opt             = True
opt@debug       = False

; We know that the shape names are stored in the array NAME
opt@shape_var   = "NAME"
opt@shape_names = (/"Colorado"/)
;opt@shape_names = (/"Colorado","Utah"/) ; You can plot more than one state if you want!

; Here we will read t2max and set all data values to missing except for
;   those over our region of interest
var_mask = shapefile_mask_data(t2max,shp_filename,opt)

; Here we will set our plot options.  Our plan is to actually make two plots on
;   top of each other.  One will plot that data for the shapes we care about, and
;   the other will plot the boundaries of all the shapes in the file
pltres = True
pltres@PanelPlot = True  ; Makes it so we can make two plots on top of each other

; Set the options for the map portion of the plot
mpres = True
; Here we can zoom in on the part of the US we are interested in.  If we don't, it will
;   just plot the whole model domain, which is quite large
; We are finding the i and j values that span the area we care about
mpres@ZoomIn = True
var1D = ndtooned(var_mask)
IND = ind_resolve( ind(var1D.ne.var_mask@_FillValue), dimsizes(var_mask))
mpres@ZoomIn = True
x_start = min(IND(:,1)) - 5
y_start = min(IND(:,0)) - 5
y_end = max(IND(:,0)) + 5
x_end = max(IND(:,1)) + 5
mpres@Xstart = x_start
mpres@Ystart = y_start
mpres@Xend = x_end
mpres@Yend = y_end
delete (var1D)
delete(IND)

; Here we apply those i and j values to the t2max data
var_mask_zoom = var_mask(y_start:y_end,x_start:x_end)

; Options for our data contour
opts = True
opts@cnFillOn = True  ; Setting the contour color fill
; wrf_contour draw contours for the t2max zoomed data
contour_mask = wrf_contour(a,wks,var_mask_zoom,opts)

; Here we are plotting the t2max contoured data over a map using the function wrf_map_overlays
plot_mask = wrf_map_overlays(a,wks,contour_mask,pltres,mpres)
; Here is where we make the second plot that draw the boundaries of all the shapes in the file
id_mask = gsn_add_shapefile_polylines(wks,plot_mask,shp_filename,True)

; Draw the plot
draw(plot_mask)
; Advance the frame
frame(wks)

; End of program
end