Skip to content

Instantly share code, notes, and snippets.

@SwampThingPaul

SwampThingPaul/Champ_cal.R

Last active August 27, 2025 16:31
Show Gist options
  • Save SwampThingPaul/09806627a0501fafbca12d410e188d3d to your computer and use it in GitHub Desktop.
Save SwampThingPaul/09806627a0501fafbca12d410e188d3d to your computer and use it in GitHub Desktop.
Download ZIP
Lake Champlain Calendar Plot
Raw
Champ_cal.R
## Title: Lake Champlain Calendar Plot
## Created by: Paul Julian ([email protected])
## Created on: 2025-08-25
# Libraries
# devtools::install_github("SwampThingPaul/AnalystHelper")
library(AnalystHelper);
library(plyr)
library(reshape2)
library(dataRetrieval)
library(zoo)
# GIS
library(sf)
library(terra)
library(tmap)
library(ggplot2)
library(ggspatial)
library(cowplot)
theme_set(theme_minimal(base_size = 16)+
theme_bw()+
theme(panel.border = element_rect("black",fill=NA,linewidth=1)))
# Helper variables
# epsg.io
wgs84 <- st_crs("EPSG:4326")
## Functions
pretty_labels <- function(cuts, digits=2){
labs <- character(length(cuts)-1)
for(i in seq_along(labs)){
if(is.infinite(cuts[i]) & !is.infinite(cuts[i+1])){
labs[i] <- paste0("< ", formatC(cuts[i+1], format="f", digits=digits))
} else if(!is.infinite(cuts[i]) & is.infinite(cuts[i+1])){
labs[i] <- paste0("> ", formatC(cuts[i], format="f", digits=digits))
} else {
labs[i] <- paste0(formatC(cuts[i], format="f", digits=digits),
" – ",
formatC(cuts[i+1], format="f", digits=digits))
}
}
labs
}
# Spatial -----------------------------------------------------------------
link <- "https://services1.arcgis.com/BkFxaEFNwHqX3tAw/arcgis/rest/services/FS_VCGI_OPENDATA_V_WATER_LKCH5K_POLY_SP_v1/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"
champ <- link|>
st_read()|>
st_transform(wgs84)
# Data --------------------------------------------------------------------
dates <- c(as.Date("1974-10-01"),as.Date(Sys.Date())-1)
site_id <- "04294500";# LAKE CHAMPLAIN AT BURLINGTON, VT
pCode <- "62614";# Lake or reservoir water surface elevation above NGVD 1929, feet
loc.dat <- readNWISsite(site_id)
loc.dat.shp <- st_as_sf(loc.dat,coords=c("dec_long_va","dec_lat_va"),crs=wgs84)
## Stage data
dat <- readNWISdv(site_id, pCode, dates[1], dates[2])
dat$Date <- date.fun(dat$Date)
# Fill empty dates
fill <- data.frame(Date=date.fun(seq(dates[1],dates[2],"1 days")))
dat <- merge(dat,fill,"Date",all.y=T)
subset(dat,is.na(X_62614_00003))
dat$WY <- WY(dat$Date,"Fed")
dat$DOWY <- hydro.day(dat$Date,"Fed")
# gap fill
dat$STG29 <- dat$X_62614_00003
dat$inter.STG29 <- na.approx(dat$STG29)
stg.cat <- c(-Inf,quantile(dat$STG29,probs=c(0.05,0.1,0.25,0.50,0.75,0.9,0.95),na.rm=T)|>as.numeric(),Inf)
dat$stg.cat <- as.factor(findInterval(dat$inter.STG29,stg.cat))
range(dat$inter.STG29)
## explore data
subset(dat,Date == date.fun(Sys.Date()-1))$inter.STG29
subset(dat,inter.STG29<=subset(dat,Date == date.fun(Sys.Date()-1))$inter.STG29)
DOWY_yest <- subset(dat,DOWY == hydro.day(date.fun(Sys.Date()-1),"Fed"))
DOWY_yest$rank.STG <- rank(DOWY_yest$inter.STG29)
DOWY_yest[order(DOWY_yest$rank.STG),]
plot(STG29~WY,DOWY_yest)
# Calendar Plot -----------------------------------------------------------
xlabs <- seq(date.fun("2025-10-01"),date.fun("2026-09-30"),"1 months");# just for months
xlabs.nums <- hydro.day(xlabs,"Fed")
xlabs.labs <- format(xlabs,"%b")
ylabs.bks <- seq(min(dat$WY),max(dat$WY),1)
ylabs <- ylabs.bks
ylabs[-seq(1,length(ylabs.bks),5)] <- " "
ylims <- range(dat$WY)+c(-0.5,0.5)
# cols <- hcl.colors(length(stg.cat)-1,palette="Spectral")
cols <- viridis::turbo(length(stg.cat)-1,direction = -1); # adjusted for more color-blind friendly color palette
champ.calendar <- ggplot(dat, aes(x = DOWY, y = WY, fill = stg.cat)) +
geom_tile(aes(group = stg.cat), colour = NA)+
scale_y_reverse(expand = c(0, 0),
limits = rev(ylims),
labels = ylabs,
breaks = ylabs.bks) +
scale_x_continuous(expand = c(0, 0),
limits = c(1,365),
breaks = xlabs.nums,
labels = xlabs.labs)+
scale_fill_manual(values = cols,
name="Lake Stage Category\n(Ft, NGVD29)",
breaks=seq(1,length(stg.cat)-1,1),
labels=pretty_labels(stg.cat, digits=2)) +
labs(title = "Lake Champlain at Burlington, VT",
subtitle = "Stage Elevation (NGVD29)",
caption = paste0("Produced: ",format(Sys.Date(),"%d %b %Y")),
x="Day of Water Year",
y="Water Year") +
theme_bw()
champ.calendar.noleg <- champ.calendar + theme(legend.position = "none")
champ.map <- ggplot() +
geom_sf(data = champ, fill = "lightblue", color = "black",linewidth =0.25) +
geom_sf(data = loc.dat.shp, color = "red", size = 3.5)+
theme_void()
legend_plot <- get_legend(champ.calendar)
champ.cal.plto <- ggdraw()+
draw_plot(champ.calendar.noleg, x = 0, y = 0, width = 0.65, height = 1)+
draw_plot(champ.map, x = 0.65, y = 0.55, width = 0.35, height = 0.45) + # inset top-right
draw_plot(legend_plot, x = 0.65, y = 0.2, width = 0.35, height = 0.25) + # legend bottom-right
theme(plot.background = element_rect(fill = "white", color = NA))
# ggsave("champ_cal_plot.png",champ.cal.plto,
# width = 6.5, height = 6, units = "in", dpi = 200)
# END ---------------------------------------------------------------------
Raw
Champ_trib.R
## Title: Lake Champlain Tributaries Plot
## Created by: Paul Julian ([email protected])
## Created on: 2025-08-27
# Libraries
# devtools::install_github("SwampThingPaul/AnalystHelper")
library(AnalystHelper);
library(plyr)
library(reshape2)
library(dataRetrieval)
library(zoo)
# GIS
library(sf)
library(terra)
# Helper variables
# epsg.io
wgs84 <- st_crs("EPSG:4326")
## Functions
pretty_labels <- function(cuts, digits=2){
labs <- character(length(cuts)-1)
for(i in seq_along(labs)){
if(is.infinite(cuts[i]) & !is.infinite(cuts[i+1])){
labs[i] <- paste0("< ", formatC(cuts[i+1], format="f", digits=digits))
} else if(!is.infinite(cuts[i]) & is.infinite(cuts[i+1])){
labs[i] <- paste0("> ", formatC(cuts[i], format="f", digits=digits))
} else {
labs[i] <- paste0(formatC(cuts[i], format="f", digits=digits),
" – ",
formatC(cuts[i+1], format="f", digits=digits))
}
}
labs
}
shade_quantiles <- function(x_range, quant_vals, cols = NULL) {
# default colors if none supplied
if (is.null(cols)) {
cols <- viridis::turbo(length(quant_vals) - 1, direction = -1)
}
# draw polygons between adjacent quantiles
for (i in seq_len(length(quant_vals) - 1)) {
polygon(
x = c(x_range[1], x_range[2], x_range[2], x_range[1]),
y = c(quant_vals[i], quant_vals[i], quant_vals[i+1], quant_vals[i+1]),
col = cols[i],
border = NA
)
}
}
# Spatial -----------------------------------------------------------------
link <- "https://services1.arcgis.com/BkFxaEFNwHqX3tAw/arcgis/rest/services/FS_VCGI_OPENDATA_V_WATER_LKCH5K_POLY_SP_v1/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"
champ <- link|>
st_read()|>
st_transform(wgs84)
maj_subbasins_link_zip <- "https://atlas.lcbp.org/WEB/Data/LCB_2013_subbasins.zip"
maj_riv_link_zip <- "https://atlas.lcbp.org/WEB/Data/LCB_major_rivers.zip"
tmp <- tempfile(fileext = ".zip")
tmpdir <- tempdir()
download.file(maj_riv_link_zip, tmp, mode = "wb")
unzip(tmp, exdir = tmpdir)
# shp_file <- list.files(tmpdir, pattern = "\\.shp$", full.names = TRUE)
shp_file <- list.files(paste0(tmpdir,"/LCB_major_rivers"), full.names = TRUE)
maj_riv <- paste0(tmpdir,"/LCB_major_rivers")|>
st_read("major_river")|>
st_transform(wgs84)
plot(st_geometry(champ))
plot(st_geometry(maj_riv),add=T)
## Data --------------------------------------------------------------------
dates <- c(as.Date("1974-10-01"),as.Date(Sys.Date())-1)
site_ids <- c("04273500","04275500","04276500","04292500","04290500","04282500")
site_labs <- c("Saranac River","Ausable River","Boquet River","Lamoille River","Winooski River","Otter Creek")
pCode <- "00060";# Discharge, cfs
loc.dat <- readNWISsite(site_ids)
loc.dat.shp <- st_as_sf(loc.dat,coords=c("dec_long_va","dec_lat_va"),crs=wgs84)
loc.dat.shp <- loc.dat.shp|>
merge(data.frame(site_no = site_ids,map_lab = site_labs),"site_no")
plot(st_geometry(champ))
plot(st_geometry(maj_riv),add=T)
plot(st_geometry(loc.dat.shp),add=T,pch=21,bg="red",cex=2)
## discharge data
dat <- readNWISdv(site_ids, pCode, dates[1], dates[2])|>
renameNWISColumns()
dat$Date <- date.fun(dat$Date)
# Fill empty dates
fill <- expand.grid(Date=date.fun(seq(dates[1],dates[2],"1 days")),
site_no = site_ids)
dat_filled <- merge(dat,fill,c("Date","site_no"),all.y=T)
subset(dat_filled,is.na(Flow))
dat_filled$WY <- WY(dat_filled$Date,"Fed")
dat_filled$DOWY <- hydro.day(dat_filled$Date,"Fed")
dat_filled$flow_acft <- cfs.to.acftd(dat_filled$Flow)
## plot --------------------------------------------------------------------
layout(matrix(c(1,2,13,7,8,
3,4,13,9,10,
5,6,14,11,12),
3,5,byrow=T),widths=c(0.25,1,1,0.25,1))
par(oma=c(2.75,3.5,1,0.5))
cols <- viridis::turbo(6,direction = -1)
for(i in seq_along(site_ids)){
par(mar=c(1,1,0.75,0));
tmp_dat <- subset(dat_filled,site_no == site_ids[i])
rng_vals <- range(tmp_dat$flow_acft,na.rm=T)
quant_vals <- quantile(tmp_dat$flow_acft,probs=c(0.05,0.1,0.25,0.50,0.75,0.9,0.95),na.rm=T)|>as.numeric()
ylim.val <- c(rng_vals[1],quant_vals[7]*1.5);ymaj <- log.scale.fun(ylim.val,"major");ymin <- log.scale.fun(ylim.val,"minor")
boxplot(tmp_dat$flow_acft,outline=F,log="y",ylim=ylim.val,ann=F,axes=F,col=NA,border=NA)
abline(h=ymaj,lty=3,col="grey")
boxplot(tmp_dat$flow_acft,outline=F,log="y",ylim=ylim.val,ann=F,axes=F,add=T)
axis_fun(2,ymaj,ymin,ymaj);box(lwd=1)
mtext(side=3,adj=0,site_labs[i],cex=0.75)
par(mar=c(1,0,0.75,0.5))
xlim.val <- date.fun(c(as.Date(Sys.Date())-90,as.Date(Sys.Date())-1));
xmaj <- seq(xlim.val[1],xlim.val[2],"30 days");
xmin <- seq(xlim.val[1],xlim.val[2],"1 days")
plot(Flow~Date,tmp_dat,ylim=ylim.val,xlim=xlim.val,ann=F,axes=F,type="n",log="y")
abline(h=ymaj,v=xmaj,lty=3,col="grey")
xx <- date.fun(c(as.Date(Sys.Date())-94,as.Date(Sys.Date())+3))
shade_quantiles(xx,quant_vals,cols = adjustcolor(cols,0.5))
lines(flow_acft~Date,tmp_dat,col="dodgerblue1",lwd=2)
#axis_fun(2,ymaj,ymin,ymaj);
box(lwd=1)
if(i%in%c(3,6)){axis_fun(1,xmaj,xmin,format(xmaj,"%m-%d"),line=-0.5)}else{axis_fun(1,xmaj,xmin,NA)}
if(i%in%c(3,6)){mtext(side=1,line=2,"Date (Month-Day)")}
}
mtext(side=2,outer=T,line=1.5,expression("Discharge (Ac-Ft d"^" -1"*")"))
bbox.lims <- st_bbox(champ)
par(mar=c(1,0,0.75,2))
plot(st_geometry(champ),ylim=bbox.lims[c(2,4)],xlim=bbox.lims[c(1,3)],
col="lightblue",border="grey",lwd=0.5,bg=NA)
plot(st_geometry(maj_riv),add=T,col="lightblue")
plot(st_geometry(loc.dat.shp),add=T,pch=21,bg="red",cex=1.5,lwd=0.5)
st_txt(loc.dat.shp[c(1,3),],loc.dat.shp$map_lab[c(1,3)],cex=0.75,pos=2)
st_txt(loc.dat.shp[2,],loc.dat.shp$map_lab[2],cex=0.75,pos=3)
st_txt(loc.dat.shp[c(4:6),],loc.dat.shp$map_lab[c(4:6)],cex=0.75,pos=4)
mapmisc::scaleBar(champ,"bottom",bty="n",cex=1,outer=F)
plot(0:1,0:1,ann=F,axes=F,type="n")
title.val <- paste0("Period of Record Quantile Ranges\n(",format(dates[1],"%d %b %Y")," - ",format(dates[2],"%d %b %Y"),")")
legend("center",legend = pretty_labels(c(0.05,0.1,0.25,0.50,0.75,0.9,0.95)),
pch=c(22),pt.bg= adjustcolor(cols,0.5),
lty=c(NA),lwd=c(0.1),col=NA,
pt.cex=2,cex=1,ncol=1,bty="n",y.intersp=1,x.intersp=0.75,xpd=NA,xjust=0.5,yjust=0.5,
title = title.val)
mtext(side=1,paste0("Produced: ",format(Sys.Date(),"%d %b %Y")),cex=0.75)
# END ---------------------------------------------------------------------
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment