HRRR Model (cutout) Domain
hname <- "hrrr201834418.nc"
#tidync(hname)
hlat <- read_ncdf(hname, var = "latitude")
hlon <- read_ncdf(hname, var = "longitude")-360
ll <- setNames(c(hlon, hlat), c("x", "y"))
st_htem <- read_ncdf(hname, var = "APCP_surface") %>%
st_as_stars(curvilinear = ll)
st_htem.proj <- st_transform(st_htem, st_crs("+init=epsg:4326"))
plot(st_htem, breaks = "equal", reset = FALSE, col = brewer.pal(n = 9, name = "Blues"),
axes = TRUE, as_points = F, border = NA, downsample = T)

Lake Champlain Watershed
library(rgdal)
rBasin <- readOGR(dsn = "LC_polygons", layer = "lake_champlain_basin_simplified", verbose = F)
rBasin.proj <- spTransform(rBasin, CRSobj = "+init=epsg:4326")
plot(rBasin, axes = TRUE, border = "red")

Add basemaps
library(maps)
#pdf("hrrr.pdf")
plot(st_htem, breaks = "equal", reset = FALSE, col = brewer.pal(n = 9, name = "Blues"),
axes = TRUE, as_points = F, border = NA, downsample = T)
#map('world', add = TRUE, col = 'black')
map('state', add = TRUE, col = 'black')
map('state', regions = c("New York", "Vermont"), add = TRUE, col = 'blue')
plot(rBasin, add = T, border = 'red')

#dev.off()