# phenocamr shiny interface / server back-end
# see matching ui.R code for formatting
# location of the data to be used in
# strings throughout
path = tempdir()
# set unix time origin
unix = "1970-01-01"
# list all colours to be used in plotting
# this allows for quick changes
col_sos_10 = "#7FFF00"
col_sos_25 = "#66CD00"
col_sos_50 = "#458B00"
col_eos_10 = "#8B6508"
col_eos_25 = "#CD950C"
col_eos_50 = "#FFB90F"
col_ci = "#C8C8C8"
col_line = "#787878"
col_text = "#787878"
# grab the latest roi list using jsonlite
# this should work across all platforms regardless
df = list_rois()
df = df[, c(
"site",
"veg_type",
"roi_id_number",
"first_date",
"last_date",
"site_years",
"lat",
"lon",
"description",
"missing_data_pct"
)]
# download metadata, and select useful columns for the explorer
metadata = list_sites()
# use gridded (daymet / worldclim) data in case of missing site specific MAT / MAP
metadata$MAT = ifelse(is.na(metadata$MAT_site),
metadata$MAT_worldclim,
metadata$MAT_site)
metadata$MAP = ifelse(is.na(metadata$MAP_site),
metadata$MAP_worldclim,
metadata$MAP_site)
# what variables do I retain
metadata = metadata[, c("site",
"ecoregion",
"koeppen_geiger",
"landcover_igbp",
"elev",
"MAT",
"MAP")]
# merge the roi list with the short metadata based on sitename
df = merge(df, metadata, by = "site")
# introduce jitter on lat/long coordinates
# this avoids that site locations with different PFTs
# will overlap, although not geographically accurate
# it helps visualize things
df$lat_j = df$lat + rnorm(length(df$lat)) * 0.00005
df$lon_j = df$lon + rnorm(length(df$lat)) * 0.00005
# subset data to exclude certain PFT / ROI classes which are irrelevant
# (no vegetation, bad ROI, mixed data types etc)
df = df[-which(
df$veg_type == "XX" | df$veg_type == "MX" |
df$veg_type == "UN" | df$veg_type == "NV" |
df$veg_type == "RF"
), ]
row.names(df) = 1:dim(df)[1]
# create a character field with html to call as a marker
# popup. This includes a thumbnail of the site.
df$preview = unlist(lapply(df$site, function(x)
paste(
"<table width=150px, border=0px>",
"<tr>",
"<td><b>",
x,
"</b></td>",
"</tr>",
"</table>",
sep = ""
)))
# start server routine
server = function(input, output, session) {
# Reactive expression for the data subsetted
# to what the user selected
v1 = reactiveValues()
v2 = reactiveValues()
reset = reactiveValues()
# the list holding the the filenames and sizes of
# the icons in the overview map
vegIcons = iconList(
WL = makeIcon("wetland_o.png", "wetland_o.png", 18, 18),
DB = makeIcon("broadleaf_o.png", "broadleaf_o.png", 18, 18),
EN = makeIcon("conifer_o.png", "conifer_o.png", 18, 18),
DN = makeIcon("conifer_o.png", "conifer_o.png", 18, 18),
EB = makeIcon("tropical_o.png", "tropical_o.png", 18, 18),
GR = makeIcon("grass_o.png", "grass_o.png", 18, 18),
SH = makeIcon("grass_o.png", "grass_o.png", 18, 18),
TN = makeIcon("grass_o.png", "grass_o.png", 18, 18),
AG = makeIcon("agriculture_o.png", "agriculture_o.png", 18, 18)
)
# function to subset the site list based upon coordinate locations
filteredData = function() {
if (!is.null(isolate(v2$lat))) {
if (input$colors == "ALL") {
tmp = df[which(
df$lat < isolate(v1$lat) &
df$lat > isolate(v2$lat) &
df$lon > isolate(v1$lon) & df$lon < isolate(v2$lon)
),]
unique(tmp)
} else{
df[which(
df$lat < isolate(v1$lat) &
df$lat > isolate(v2$lat) &
df$lon > isolate(v1$lon) &
df$lon < isolate(v2$lon) & df$veg_type == input$colors
),]
}
} else{
if (input$colors == "ALL") {
unique(df)
} else{
df[df$veg_type == input$colors,]
}
}
}
getValueData = function(table) {
# calculate the total number of site years
# and sites in the dataset
total_site_years = sum(table$site_years, na.rm = TRUE)
total_sites = length(unique(table$site))
output$site_count = renderInfoBox({
infoBox("Sites",
total_sites,
icon = icon("list"),
color = "blue",
fill = TRUE
)
})
output$year_count = renderInfoBox({
infoBox("Site Years",
total_site_years,
icon = icon("list"),
color = "blue",
fill = TRUE
)
})
}
# fill site count etc fields
getValueData(df)
# create function to plot the climatology
# occurs too many times to repeat the code
updateClimatology = function(){
output$climate = renderPlot({
par(mar = c(4, 4, 4, 1))
plot(
1,
1,
type = 'n',
xaxt = 'n',
yaxt = 'n',
xlab = '',
ylab = '',
bty = 'n',
xlim = c(0, 100),
ylim = c(0, 100)
)
plot(MAP~MAT,data=filteredData(),
xlab=expression("MAT ("*degree*"C)"),
ylab="MAP (mm)",
pch=19,
col=rgb(0.5,0.5,0.5,0.3),
xlim=c(-15,30),
ylim=c(0,3000),
main = "Climatology"
)
}, height = function() {
session$clientData$output_climate_height
})
}
# Use leaflet() here, and only include aspects of the map that
# won't need to change dynamically (at least, not unless the
# entire map is being torn down and recreated).
output$map = renderLeaflet({
map = leaflet(df) %>%
addTiles(
"http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.jpg",
attribution = 'Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',
group = "World Imagery"
) %>%
addWMSTiles(
"http://webmap.ornl.gov/ogcbroker/wms?",
layers = "10004_31",
options = WMSTileOptions(format = "image/png", transparent = TRUE),
attribution = "MODIS Land Cover (MCD12Q1) © NASA",
group = "MODIS Land Cover"
) %>%
addProviderTiles(
"OpenTopoMap",
group = "Open Topo Map"
) %>%
addMarkers(
lat = ~ lat_j,
lng = ~ lon_j,
icon = ~ vegIcons[veg_type],
popup = ~ preview
) %>%
# Layers control
addLayersControl(
baseGroups = c("World Imagery","MODIS Land Cover","Open Topo Map"),
position = c("topleft"),
options = layersControlOptions(collapsed = TRUE)
) %>%
setView(lng = 11,
lat = 45,
zoom = 2)
})
# Incremental changes to the map. Each independent set of things that can change
# should be managed in its own observer.
observe({
leafletProxy("map", data = filteredData()) %>%
clearMarkers() %>%
addMarkers(
lat = ~ lat_j,
lng = ~ lon_j,
icon = ~ vegIcons[veg_type],
popup = ~ preview
)
# update the data table in the explorer
output$table = DT::renderDataTable({
tmp = filteredData()#[,-DT_drop]
tmp = tmp[,-c((ncol(tmp)-2):ncol(tmp))]
return(tmp)
},
selection = "single",
options = list(lengthMenu = list(c(5, 10), c('5', '10')),
pom = list('site_years')),
extensions = c('Responsive'))
# update value box
getValueData(filteredData())
# update climatology plot
updateClimatology()
})
# grab the bounding box, by clicking the map
observeEvent(input$map_click, {
# if clicked once reset the bounding box
# and show all data
if (!is.null(isolate(v2$lat))) {
# set bounding box values to NULL
v1$lat = NULL
v2$lat = NULL
v1$lon = NULL
v2$lon = NULL
leafletProxy("map", data = filteredData()) %>%
clearMarkers() %>%
clearShapes() %>%
addMarkers(
lat = ~ lat_j,
lng = ~ lon_j,
icon = ~ vegIcons[veg_type],
popup = ~ preview
)
getValueData(filteredData())
# update climatology plot
updateClimatology()
} else{
# grab bounding box coordinates
# TODO: validate the topleft / bottom right order
if (!is.null(isolate(v1$lat))) {
v2$lat = input$map_click$lat
v2$lon = input$map_click$lng
} else{
v1$lat = input$map_click$lat
v1$lon = input$map_click$lng
leafletProxy("map", data = filteredData()) %>%
clearMarkers() %>%
addMarkers(
lat = ~ lat_j,
lng = ~ lon_j,
icon = ~ vegIcons[veg_type],
popup = ~ preview
) %>%
addCircleMarkers(
lng = isolate(v1$lon),
lat = isolate(v1$lat),
color = "red",
radius = 3,
fillOpacity = 1,
stroke = FALSE
)
}
}
# if the bottom right does exist
if (!is.null(isolate(v2$lat))) {
# subset data based upon topleft / bottomright
# first put all data in tmp data table
tmp = filteredData()
# check if the dataset is not empty
if (dim(tmp)[1] != 0) {
# update the map
leafletProxy("map", data = tmp) %>%
clearMarkers() %>%
addMarkers(
lat = ~ lat_j,
lng = ~ lon_j,
icon = ~ vegIcons[veg_type],
popup = ~ preview
) %>%
addRectangles(
lng1 = isolate(v1$lon),
lat1 = isolate(v1$lat),
lng2 = isolate(v2$lon),
lat2 = isolate(v2$lat),
fillColor = "transparent",
color = "red"
)
# update the data table in the explorer
output$table = DT::renderDataTable({
tmp = filteredData()#[,-DT_drop]
tmp = tmp[,-c((ncol(tmp)-2):ncol(tmp))]
return(tmp)
},
selection = "single",
options = list(lengthMenu = list(c(5, 10), c('5', '10')),
pom = list('site_years')),
extensions = c('Responsive'))
# update the value box
getValueData(filteredData())
# update climatology plot
updateClimatology()
} else{
# set bounding box values to NULL
v1$lat = NULL
v2$lat = NULL
v1$lon = NULL
v2$lon = NULL
leafletProxy("map", data = filteredData()) %>%
clearMarkers() %>%
clearShapes() %>%
addMarkers(
lat = ~ lat_j,
lng = ~ lon_j,
icon = ~ vegIcons[veg_type],
popup = ~ preview
)
# then stuff it back into the shiny output
# object
output$table = DT::renderDataTable({
tmp = filteredData()
tmp = tmp[,-c((ncol(tmp)-2):ncol(tmp))]
return(tmp)
},
selection = "single",
options = list(lengthMenu = list(c(5, 10), c('5', '10')),
pom = list('site_years')),
extensions = c('Responsive'))
# update climatology plot
updateClimatology()
}
}
})
downloadData = function(myrow, frequency, percentile) {
# if nothing selected return NULL
if (length(myrow) == 0) {
return(NULL)
}
# Create a Progress object
# Make sure it closes when we exit this reactive, even if there's an error
progress = shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Making plot", value = 0)
# grab the necessary parameters to download the site data
site = as.character(df$site[as.numeric(myrow)])
veg_type = as.character(df$veg_type[as.numeric(myrow)])
roi_id = as.numeric(df$roi_id_number[as.numeric(myrow)])
# formulate the filename of the data we need to download
data_file = sprintf("%s/%s_%s_%04d_%sday.csv",
path,
site,
veg_type,
roi_id,
frequency)
# check if the data is on file, if so assume it is
# checked for outliers and has smoothed data
if (file.exists(data_file)) {
# read the data to plot_data
data_raw = read_phenocam(data_file)
data = data_raw$data
} else{
# kick start progress bar
progress$set(value = 0.5, detail = "Downloading PhenoCam data")
# download phenocam data from the server
# do not detect outliers and smooth in this step
# do this separately to allow for proper progresss
# bar updating (keep people busy as this takes a while)
download_phenocam(
site = site,
veg_type = veg_type,
roi_id = roi_id,
frequency = input$frequency,
outlier_detection = TRUE,
smooth = TRUE,
out_dir = path
)
# read the data to plot_data
data_raw = read_phenocam(data_file)
data = data_raw$data
}
# grab the transition dates
# trap errors, mainly if no dates can be detected return an empty
# string (NAs) to prevent plotting errors further down.
# code this up in the transition.dates() function TODO TODO TODO
spring = transition_dates(data_raw,
percentile = percentile,
reverse = FALSE)
fall = transition_dates(data_raw,
percentile = percentile,
reverse = TRUE)
# Final plot preparations
progress$set(value = 0.8, detail = "preparing final plot")
# formulate variable names dynamically
gcc_val = sprintf("gcc_%s", percentile)
rcc_val = sprintf("rcc_%s", percentile)
#out_val = data[,sprintf("outlierflag_gcc_%s", percentile)]
col_val = colnames(data)
smooth_val_gcc = sprintf("smooth_gcc_%s", percentile)
smooth_val_rcc = sprintf("smooth_rcc_%s", percentile)
ci_val = sprintf("smooth_ci_gcc_%s", percentile)
# stuff things in reactive value
date = data$date
doy = data$doy
year = data$year
# calculate gcc / bcc / rcc time series
gcc = data[, which(col_val == gcc_val)]
bcc = data$midday_b / (data$midday_r + data$midday_g + data$midday_b)
rcc = data$midday_r / (data$midday_r + data$midday_g + data$midday_b)
# put outlier and interpolation flags into readable
# variable names
out = data[,sprintf("outlierflag_gcc_%s", percentile)]
int_flag = data$int_flag
# rename outlier values with proper descriptive names
gcc_smooth = data[, which(col_val == smooth_val_gcc)]
ci = data[, which(col_val == ci_val)]
# create data frame which is exported and used in
# all plotting routines
plot_data = data.frame(date, year, doy, gcc, rcc, bcc, out, gcc_smooth, ci, int_flag)
# return raw data, and derived phenology metric
# in a structured list (split up later on depending on use)
return(list(plot_data, spring, fall))
}
# observe the state of the table, if changed update the data
inputData = reactive({
downloadData(
as.numeric(input$table_row_last_clicked),
as.numeric(input$frequency),
input$percentile
)
})
# plot the data ---------------
output$time_series_plot = renderPlotly({
# grab plotting data
data = inputData()
# grab some site specific data (sitename etc)
myrow = as.numeric(input$table_row_last_clicked)
site = as.character(df$site[as.numeric(myrow)])
veg = as.character(df$veg_type[as.numeric(myrow)])
roi = as.numeric(df$roi_id_number[as.numeric(myrow)])
# axis styling
ax = list(
title = "",
zeroline = FALSE,
showline = FALSE,
showticklabels = FALSE,
showgrid = FALSE
)
# plotting routine
if (is.null(data) | all((is.na(data[[1]])) & all(is.na(data[[2]]))) ) {
# populate the download handler with error message
output$downloadData <- downloadHandler(
filename = sprintf(
'phenology_data_%s_%s_%04d_%s.csv',
site,
veg,
roi,
Sys.Date()
),
content = function(file = filename) {
write.table(
"NO DATA",
file,
row.names = FALSE,
col.names = TRUE,
quote = FALSE,
sep = ','
)
}
)
# plot error message
p = plot_ly(
x = 0,
y = 0,
text = "NO DATA - SELECT A (DIFFERENT) SITE",
mode = "text",
textfont = list(color = '#000000', size = 16)
) %>%
layout(xaxis = ax, yaxis = ax)
} else{
# split up the data into true plotting data (time series) and annotation (phenology dates)
plot_data = data[[1]]
spring = data[[2]]
fall = data[[3]]
# format dates correctly (unix time to date format)
spring[, 1:9] = apply(spring[, 1:9], 2, function(x)
as.character(as.Date(x, origin = unix)))
fall[, 1:9] = apply(fall[, 1:9], 2, function(x)
as.character(as.Date(x, origin = unix)))
# bind spring and fall phenology data in a coherent format
phenology = rbind(spring,fall)
rising_length = dim(spring)[1]
falling_length = dim(fall)[1]
direction = c(rep("rising", rising_length),
rep("falling", falling_length))
sitename = rep(site,rising_length + falling_length)
veg_type = rep(veg,rising_length + falling_length)
roi_id = rep(sprintf("%04d",roi),rising_length + falling_length)
phenology[, 1:9] = apply(phenology[, 1:9], 2, function(x)
as.character(as.Date(x, origin = unix)))
# bind in new labels
phenology = cbind(sitename,veg_type,roi_id,direction,phenology)
# drop NA lines
phenology = na.omit(phenology)
# create header with output information
phenology_header = matrix("#",12,1)
# populate the header file
phenology_header[2,1] = "# Transition date estimates"
phenology_header[4,1] = sprintf("# Sitename: %s",site)
phenology_header[5,1] = sprintf("# Vegetation Type: %s",veg)
phenology_header[6,1] = sprintf("# ROI ID: %04d",roi)
phenology_header[7,1] = sprintf("# Aggregation period: %s", as.numeric(input$frequency))
phenology_header[8,1] = sprintf("# Year min: %s",
min(strptime(as.matrix(phenology[, 5:13]),"%Y-%m-%d")$year +
1900), # THIS IS UGLY THIS CAN BE SHORTER
na.rm = TRUE)
phenology_header[9,1] = sprintf("# Year max: %s",
max(strptime(as.matrix(phenology[, 5:13]),"%Y-%m-%d")$year +
1900),
na.rm = TRUE)
phenology_header[10,1] = sprintf("# Creation Date: %s", Sys.Date())
phenology_header[11,1] = sprintf("# Creation Time: %s", format(Sys.time(), "%H:%M:%S"))
# populate the download handler with phenology data
output$downloadData <- downloadHandler(
filename = function(){sprintf(
'phenology_data_%s_%s_%04d_%s.csv',
site,
veg,
as.numeric(roi),
Sys.Date()
)},
content = function(file = filename) {
write.table(
phenology_header,
file,
row.names = FALSE,
col.names = FALSE,
quote = FALSE
)
write.table(
phenology,
file,
row.names = FALSE,
col.names = TRUE,
quote = FALSE,
sep = ',',
append = TRUE
)
}
)
# add an ID column to the phenology data
# for easy reshaping and plotting
spring = data.frame(id=1:nrow(spring),spring)
fall = data.frame(id=1:nrow(fall),fall)
# full time series
if (input$plot_type == "bydate") {
# duplicate outlier flag data, for plotting
plot_data$outlier_symbols = plot_data$out
plot_data$colours = plot_data$out
plot_data$colours[plot_data$colours == 0] = rgb(0.3, 0.3, 0.3)
plot_data$colours[plot_data$colours == 1] = rgb(1, 0, 0)
# rename the outliers symbols
plot_data$outlier_symbols[plot_data$outlier_symbol == 0] = "circle"
plot_data$outlier_symbols[plot_data$outlier_symbol == 1] = "circle-open"
# convert date
plot_data$date = as.Date(plot_data$date,"%Y-%m-%d")
# when showing multiple time series
# remove the transition dates
if (input$rccbcc == "FALSE") {
# base plot shows the Gcc confidence interval
# switch order for cleaner plotting
p = plot_ly(
data = plot_data,
x = ~ date,
y = ~ gcc_smooth,
showlegend = FALSE
) %>%
add_trace(
x = ~ date,
y = ~ gcc_smooth - ci,
mode = "lines",
type = "scatter",
line = list(width = 0),
showlegend = FALSE,
name = "Gcc 95% CI"
) %>%
add_trace(
y = ~ gcc_smooth + ci,
fill = "tonexty",
mode = "lines",
line = list(width = 0),
color = I(rgb(0.5,0.5,0.5)),
showlegend = TRUE,
name = "Gcc 95% CI"
) %>%
add_trace(
y = ~ gcc_smooth,
mode = "lines",
line = list(width = 2, color = rgb(0.3,0.3,0.3)),
name = "Gcc loess fit",
showlegend = TRUE
) %>%
add_trace(
y = ~ gcc,
mode = "markers",
type = "scatter",
symbol = ~ I(outlier_symbols),
color = ~ I(rep(rgb(0.3, 0.3, 0.3),nrow(plot_data))),
name = "Gcc",
showlegend = TRUE
) %>%
# SOS spring
# 10%
add_trace(
data = spring,
x = ~ as.Date(transition_10),
y = ~ threshold_10,
mode = "markers",
type = "scatter",
marker = list(color = "#7FFF00", symbol = "circle"),
name = "SOS (10%)",
showlegend = TRUE
) %>%
add_segments(x = ~ as.Date(transition_10_lower_ci),
xend = ~ as.Date(transition_10_upper_ci),
y = ~ threshold_10,
yend = ~ threshold_10,
line = list(color = "#7FFF00"),
name = "SOS (10%) - CI"
) %>%
# 25 %
add_trace(
x = ~ as.Date(transition_25),
y = ~ threshold_25,
mode = "markers",
type = "scatter",
marker = list(color = "#66CD00", symbol = "square"),
showlegend = TRUE,
name = "SOS (25%)"
) %>%
add_segments(x = ~ as.Date(transition_25_lower_ci),
xend = ~ as.Date(transition_25_upper_ci),
y = ~ threshold_25,
yend = ~ threshold_25,
line = list(color = "#66CD00"),
name = "SOS (25%) - CI"
) %>%
# 50 %
add_trace(
x = ~ as.Date(transition_50),
y = ~ threshold_50,
mode = "markers",
type = "scatter",
marker = list(color = "#458B00", symbol = "diamond"),
showlegend = TRUE,
name = "SOS (50%)"
) %>%
add_segments(x = ~ as.Date(transition_50_lower_ci),
xend = ~ as.Date(transition_50_upper_ci),
y = ~ threshold_50,
yend = ~ threshold_50,
line = list(color = "#458B00"),
name = "SOS (50%) - CI"
) %>%
# EOS fall
# 50%
add_trace(
data = fall,
x = ~ as.Date(transition_50),
y = ~ threshold_50,
mode = "markers",
type = "scatter",
marker = list(color = "#FFB90F", symbol = "diamond"),
showlegend = TRUE,
name = "EOS (50%)"
) %>%
add_segments(x = ~ as.Date(transition_50_lower_ci),
xend = ~ as.Date(transition_50_upper_ci),
y = ~ threshold_50,
yend = ~ threshold_50,
line = list(color = "#FFB90F"),
name = "EOS (50%) - CI"
) %>%
# 25 %
add_trace(
x = ~ as.Date(transition_25),
y = ~ threshold_25,
mode = "markers",
type = "scatter",
marker = list(color = "#CD950C", symbol = "square"),
showlegend = TRUE,
name = "EOS (25%)"
) %>%
add_segments(x = ~ as.Date(transition_25_lower_ci),
xend = ~ as.Date(transition_25_upper_ci),
y = ~ threshold_25,
yend = ~ threshold_25,
line = list(color = "#CD950C"),
name = "EOS (25%) - CI"
) %>%
# 10 %
add_trace(
x = ~ as.Date(transition_10),
y = ~ threshold_10,
mode = "markers",
marker = list(color = "#8B6508", symbol = "circle"),
showlegend = TRUE,
name = "EOS (10%)"
) %>%
add_segments(x = ~ as.Date(transition_10_lower_ci),
xend = ~ as.Date(transition_10_upper_ci),
y = ~ threshold_10,
yend = ~ threshold_10,
line = list(color = "#8B6508"),
name = "EOS (10%) - CI"
) %>%
layout(xaxis = list(title = "Date"),
yaxis = list(title = "Gcc"))
} else {
p = plot_ly(
data = plot_data,
x = ~ date,
y = ~ bcc,
showlegend = FALSE
) %>%
add_trace(
y = ~ bcc,
mode = "markers",
type = "scatter",
marker = list(color = "#0000FF", symbol = "circle"),
name = "Bcc",
showlegend = TRUE
) %>%
add_trace(
y = ~ rcc,
marker = list(color = "#FF0000", symbol = "circle"),
name = "Rcc",
showlegend = TRUE
) %>%
add_trace(
y = ~ gcc,
symbol = ~ I(outlier_symbols),
marker = list(color = "#00FF00", symbols = ~ outlier_symbols),
name = "Gcc",
showlegend = TRUE
) %>%
add_trace(
y = ~ gcc_smooth,
type="scatter",
mode = "lines",
line = list(width = 2, color = "#00FF00"),
name = "Gcc loess fit",
showlegend = TRUE
) %>%
layout(xaxis = list(title = "Date"),
yaxis = list(title = "Gcc / Rcc / Bcc"))
}
} else{
# condense to one plot along DOY
if (input$plot_type == "byyear") {
ltm = plot_data %>% group_by(doy) %>%
summarise(mn = mean(gcc_smooth), sd = sd(gcc_smooth), doymn = mean(doy))
# delete interpolated data
# NA values == data available
plot_data$gcc_smooth[!is.na(plot_data$int_flag)] = NA
p = plot_ly(
data = ltm,
x = ~ doymn,
y = ~ mn - sd,
mode = "lines",
fill = "none",
type = 'scatter',
line = list(width = 0),
showlegend = FALSE,
name = "1 SD"
) %>%
add_trace(
y = ~ mn + sd,
mode = "lines",
fill = "tonexty",
line = list(width = 0),
color = I(rgb(0.5,0.5,0.5)),
showlegend = TRUE,
name = "1 SD"
) %>%
add_trace(
y = ~ mn,
line = list(
width = 2,
dash = "dot",
color = "black"
),
name = "LTM",
showlegend = TRUE
) %>%
add_trace(
data = plot_data,
x = ~ doy,
y = ~ gcc_smooth,
split = ~ year,
type = "scatter",
mode = "lines",
line = list(width = 2, color = "Set1"),
showlegend = TRUE
) %>%
layout(xaxis = list(title = "DOY"),
yaxis = list(title = "Gcc"))
} else {
# these are the transition date regression plots
if (dim(fall)[1] < 9 & dim(spring)[1] < 9) {
p = plot_ly(
x = 0,
y = 0,
text = "TOO FEW (<9) DATES FOR MEANINGFUL REGRESSION ANALYSIS",
mode = "text",
textfont = list(color = '#000000', size = 16)
) %>% layout(xaxis = ax, yaxis = ax)
} else {
# grab dates from the fall and spring matrices
fall_date = unique(as.Date(fall[,grep(sprintf("^transition_%s$",upper.thresh*100),names(fall))]))
spring_date = unique(as.Date(spring[,grep(sprintf("^transition_%s$",upper.thresh*100),names(spring))]))
fall_doy = as.numeric(format(fall_date, "%j"))
fall_year = as.numeric(format(fall_date, "%Y"))
spring_doy = as.numeric(format(spring_date, "%j"))
spring_year = as.numeric(format(spring_date, "%Y"))
# # regression stats
reg_spring = lm(spring_doy ~ spring_year)
reg_fall = lm(fall_doy ~ fall_year)
# summaries
reg_spring_sum = summary(reg_spring)
reg_fall_sum = summary(reg_fall)
# r-squared and slope
r2_spring = round(reg_spring_sum$r.squared, 2)
slp_spring = round(reg_spring_sum$coefficients[2, 1], 2)
r2_fall = round(reg_fall_sum$r.squared, 2)
slp_fall = round(reg_fall_sum$coefficients[2, 1], 2)
p1 = plot_ly(
x = spring_year,
y = spring_doy,
yaxis = "y1",
title = "PhenoCam Phenology (DOY)"
) %>%
add_trace(
x = spring_year,
y = spring_doy,
marker = list(color = "#66CD00", symbol = "square"),
mode = "markers",
type = "scatter",
name = "Spring",
yaxis = "y1"
) %>%
add_trace(type = "scatter",
x = fall_year,
y = fall_doy,
marker = list(color = "#CD950C", symbol = "square"),
mode = "markers",
type = "scatter",
name = "Autumn",
yaxis = "y1"
) %>%
add_trace(
x = spring_year,
y = reg_spring$fitted.values,
mode = "lines",
type = "scatter",
name = sprintf("R2: %s| slope: %s", r2_spring, slp_spring),
line = list(width = 2, color = "#66CD00"),
yaxis = "y1"
) %>%
add_trace(
x = fall_year,
y = reg_fall$fitted.values,
mode = "lines",
type = "scatter",
name = sprintf("R2: %s| slope: %s", r2_fall, slp_fall),
line = list(width = 2, color = "#CD950C"),
yaxis = "y1"
) %>%
layout(xaxis = list(title = "Year"),
yaxis = list(title = "DOY"))
}
}
}
}
}) # plotly action end
} # server function end
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.