Make a multi-series highchart plot of HOBO temperature logger data collected at an impaired urban stream. The monitoring site is the outfall of a large channel impoundment. Display mean daily water temperature by year, as well as daily min/max and inner-quartile range (IQR) for the period of record (POR). Add bio-criteria temperature threshold.
##load libraries
require(tidyverse); require(highcharter); require(stringr); require(lubridate); require(forcats)##set working directory and read data; one table per field season
glimpse(df2007)## Observations: 15,456
## Variables: 3
## $ Date_Time <fctr> 04/26/2007 0:00, 04/26/2007 0:15, 04/26/2007 0:30...
## $ StnRVT_CMB1 <dbl> 66.814, 66.686, 66.600, 66.472, 66.772, 66.814, 66...
## $ StnRVT_CMB2 <dbl> 66.857, 66.729, 66.643, 66.515, 66.729, 66.857, 66...##bind tables by row to new tbl_df
df_temp<-bind_rows(list(df2006,df2007,df2008, df2015))
##rename variables and remove known outlier
df_temp %>% rename(inflowTemp=StnRVT_CMB1, outfallTemp=StnRVT_CMB2) %>%
filter(outfallTemp < 95) -> df_temp_filter
glimpse(df_temp_filter)## Observations: 63,907
## Variables: 3
## $ Date_Time <chr> "04/15/2006 0:00", "04/15/2006 0:15", "04/15/2006 ...
## $ inflowTemp <dbl> 49.870, 49.781, 49.693, 49.604, 49.560, 49.516, 49...
## $ outfallTemp <dbl> 53.035, 52.947, 52.904, 52.904, 52.772, 52.772, 52...##format date and extract year and Julian day; select columns of interest
df_temp_filter %>% mutate(Date=str_sub(Date_Time, 1,10)) %>%
mutate(Date=as.Date(Date, format = "%m/%d/%Y")) %>%
mutate(year=year(Date)) %>%
mutate(yday = yday(Date))%>%
select(inflowTemp:yday) -> df_temp_filter2
head(df_temp_filter2, n=5)## # A tibble: 5 × 5
## inflowTemp outfallTemp Date year yday
## <dbl> <dbl> <date> <dbl> <dbl>
## 1 49.870 53.035 2006-04-15 2006 105
## 2 49.781 52.947 2006-04-15 2006 105
## 3 49.693 52.904 2006-04-15 2006 105
## 4 49.604 52.904 2006-04-15 2006 105
## 5 49.560 52.772 2006-04-15 2006 105##account for 2008 leap year; correct Julian day
df_temp_filter2 %>% mutate(yday=ifelse(year==2008, yday-1, yday )) -> df_temp_filter2
##filter dates to better align PORs across sample years
df_temp_filter2 %>%
filter(yday >= yday(as.Date("04-27-2006", format="%m-%d-%Y")) &
yday<=yday(as.Date("10-11-2006", format="%m-%d-%Y"))) -> df_temp_filter2##calculate Julian day min/max and IQR; check daily sample count using n()
by_yday<- group_by(df_temp_filter2, yday)
yday_stats <- summarise(by_yday,
count=n(),
p25=quantile(outfallTemp, probs=0.25),
p75=quantile(outfallTemp, probs=0.75),
min_yday=min(outfallTemp, na.rm=TRUE),
max_yday=max(outfallTemp, na.rm=TRUE))
head(yday_stats, n=5)## # A tibble: 5 × 6
## yday count p25 p75 min_yday max_yday
## <dbl> <int> <dbl> <dbl> <dbl> <dbl>
## 1 117 192 48.27200 50.47725 47.736 54.082
## 2 118 192 48.43875 51.50100 48.004 53.646
## 3 119 192 49.82600 51.07200 48.049 55.342
## 4 120 192 49.42700 52.93675 48.004 57.505
## 5 121 192 50.69800 53.98375 48.361 57.160##drop count field and join Julian day stats to timeseries
by_yday_stats<-select(yday_stats, -count)
##for now only display 2006-2008 data
by_date<- group_by(df_temp_filter2, Date)
by_date_stats <- summarise(by_date,
count = n())
df_fnl<-filter(by_date_stats, Date < "2015-01-01")
df_fnl<-mutate(df_fnl, yday=yday(Date))
##join min/max and IQR stats to timeseries
df_fnl2<-left_join(df_fnl, by_yday_stats, by="yday")
df_fnl2<-select(df_fnl2, -count)
glimpse(df_fnl2)## Observations: 480
## Variables: 6
## $ Date <date> 2006-04-27, 2006-04-28, 2006-04-29, 2006-04-30, 2006...
## $ yday <dbl> 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127...
## $ p25 <dbl> 48.27200, 48.43875, 49.82600, 49.42700, 50.69800, 51....
## $ p75 <dbl> 50.47725, 51.50100, 51.07200, 52.93675, 53.98375, 53....
## $ min_yday <dbl> 47.736, 48.004, 48.049, 48.004, 48.361, 50.268, 50.40...
## $ max_yday <dbl> 54.082, 53.646, 55.342, 57.505, 57.160, 56.036, 55.99...##rename variables and gather
dtempgather<-df_fnl2 %>%
select(Date, p25, p75, max_yday, min_yday) %>%
rename(iqr_min=p25,
iqr_max=p75,
rec_max=max_yday,
rec_min=min_yday) %>%
gather(key, value, -Date) ##wide to long format
##spread
dtempspread <- dtempgather %>%
separate(key, c("serie", "type"), sep = "_") %>%
spread(type, value) ##back to long format
head(dtempspread, n=5)## # A tibble: 5 × 4
## Date serie max min
## <date> <chr> <dbl> <dbl>
## 1 2006-04-27 iqr 50.47725 48.27200
## 2 2006-04-27 rec 54.08200 47.73600
## 3 2006-04-28 iqr 51.50100 48.43875
## 4 2006-04-28 rec 53.64600 48.00400
## 5 2006-04-29 iqr 51.07200 49.82600##format as factors
temps <- dtempspread %>%
mutate(serie = factor(serie, levels = c("rec", "iqr")),
serie = fct_recode(serie, Record = "rec", IQR = "iqr"))
##add timestamp for plotting
temps <- mutate(temps, dt = datetime_to_timestamp(Date))
head(temps, n=5)## # A tibble: 5 × 5
## Date serie max min dt
## <date> <fctr> <dbl> <dbl> <dbl>
## 1 2006-04-27 IQR 50.47725 48.27200 1.146096e+12
## 2 2006-04-27 Record 54.08200 47.73600 1.146096e+12
## 3 2006-04-28 IQR 51.50100 48.43875 1.146182e+12
## 4 2006-04-28 Record 53.64600 48.00400 1.146182e+12
## 5 2006-04-29 IQR 51.07200 49.82600 1.146269e+12##only need one Julian year for displaying IQR and min/max bands
temps<-filter(temps, Date < "2007-01-01")
temps<-select(temps, dt, serie, max, min)
temps<-mutate(temps, max=as.integer(max)) ##for plotting
temps<-mutate(temps, min=as.integer(min)) ##for plotting##for now, just display years 2006-2008
by_date<- filter(df_temp_filter2, Date < "2015-01-01")
##calculate mean temp
by_date<- group_by(by_date, Date)
by_date_mean <- summarise(by_date,
mean = as.integer(mean(outfallTemp, na.rm = TRUE)))
##plotting hack; treat date as Julian day and assign year 2006 to overplot annual series on same date axis
by_date_mean <- mutate(by_date_mean, year = as.character(year(Date)))
by_date_mean<-mutate(by_date_mean, plot_date=as.Date(paste("2006-",format(Date, "%m-%d"), sep=""), format="%Y-%m-%d"))
by_date_mean<- mutate(by_date_mean, dt = datetime_to_timestamp(plot_date))
head(by_date_mean, n=5)## # A tibble: 5 × 5
## Date mean year plot_date dt
## <date> <int> <chr> <date> <dbl>
## 1 2006-04-27 50 2006 2006-04-27 1.146096e+12
## 2 2006-04-28 50 2006 2006-04-28 1.146182e+12
## 3 2006-04-29 51 2006 2006-04-29 1.146269e+12
## 4 2006-04-30 53 2006 2006-04-30 1.146355e+12
## 5 2006-05-01 54 2006 2006-05-01 1.146442e+12## setup chart
hc <- highchart(width=650) %>%
hc_title(text = "StnRVT_CMB2 Water Temperature") %>%
hc_subtitle(text = "Maximum bio-criteria threshold displayed") %>%
hc_xAxis(type = "datetime", showLastLabel = TRUE,
dateTimeLabelFormats = list(day = "%d", month="%b"))%>%
hc_tooltip(shared = TRUE, useHTML = TRUE,
headerFormat = as.character(tags$small("{point.x: %b %d}", tags$br()))) %>%
hc_plotOptions(series = list(borderWidth = 0,groupPadding= 2, pointWidth = 5)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_yAxis(title = list(text = "Degrees Fahrenheit"),
plotLines = list(list(color = "darkgray", alpha=0.2, width = 2, value = 75.2, ##add Brook trout maximum tolerance
label=list(text = "max tolerance", style = list(color = 'black', fontWeight = 'bold')))))
## add temperature data to highchart
hc <- hc %>%
hc_add_series_df(temps, type = "columnrange", x=dt, low=min, high=max, group=serie)
##set colors
x<-hex_to_rgba(c("#ecf0f1", "#DCDCDC"), alpha=1)
hc<-hc_colors(hc, colors=c(x,"rgba(22,160,133,0.8)", "rgba(231,76,60,0.8)","rgba(41,128,185,0.8"))
##add annual mean temperature series
hc <- hc %>%
hc_add_series_df(by_date_mean, type = "line", x=dt, y=mean, group=year)