# Generates a figure showing the results of the different k600 models
library(LakeMetabolizer)
data.path = system.file('extdata', package="LakeMetabolizer")
sp.data = load.all.data('sparkling', data.path)
ts.data = sp.data$data #pull out just the timeseries data
#calculate U10 and add it back onto the original
u10 = wind.scale(ts.data)
ts.data = rmv.vars(ts.data, 'wnd', ignore.offset=TRUE) #drop old wind speed column
ts.data = merge(ts.data, u10) #merge new u10 into big dataset
#Calculate k600 using the k.cole wind-based model
k600_cole = k.cole(ts.data)
k600_crusius = k.crusius(ts.data)
ha2m2 <- 10000
kd = sp.data$metadata$averagekd
wnd.z = 10 #because we converted to u10
atm.press = 1018
lat = sp.data$metadata$latitude
lake.area = sp.data$metadata$lakearea*ha2m2
#for k.read and k.macIntyre, we need LW_net.
#Calculate from the observations we have available.
lwnet = calc.lw.net(ts.data, lat, atm.press)
ts.data = merge(ts.data, lwnet)
k600_read = k.read(ts.data, wnd.z=wnd.z, Kd=kd, atm.press=atm.press,
lat=lat, lake.area=lake.area)
k600_soloviev = k.read.soloviev(ts.data, wnd.z=wnd.z, Kd=kd,
atm.press=atm.press, lat=lat, lake.area=lake.area)
k600_macIntyre = k.macIntyre(ts.data, wnd.z=wnd.z, Kd=kd, atm.press=atm.press)
k600_heiskanen = k.heiskanen(ts.data, wnd.z, kd, atm.press)
k600_vachon = k.vachon(ts.data, lake.area)
# ---- figure ----
cols <- c("#1b9e77", "#d95f02", "black", "#e7298a", "DodgerBlue", "#e6ab02", "grey50")
models <- list(
list('name'="MacIntyre", data = k600_macIntyre, col = cols[1], lty = 6, lwd = 1.7),
list('name'="Cole", data = k600_cole, col = cols[2], lty = 1, lwd = 1.2),
list('name'="Vachon", data = k600_vachon, col = cols[3], lty = 1, lwd = 1.1),
list('name'="Read", data = k600_read, col = cols[4], lty = 1, lwd = 1.2),
list('name'="Soloviev", data = k600_soloviev, col = cols[5], lty = 6, lwd = 1.7),
list('name'="Heiskanen", data = k600_heiskanen, col = cols[6], lty = 1, lwd = 1.2),
list('name'="Crusius", data = k600_crusius, col = cols[7], lty = 1, lwd = 1.2))
add_axes <- function(xlim, ylim, ylabel = pretty(ylim,10), panel.txt){
prc_x = 0.1 # position for text relative to axes
prc_y = 0.07
tick_len <- 0.15
ext_x <- c(xlim[1]-86400, pretty(xlim,3), xlim[2]+86400)
ext_y <- c(ylim[1]-10, pretty(ylim,10), ylim[2]+10)
ylab <- c("",ylabel,"")
if (is.na(ylabel[1])) ylab = NA
axis(side = 1, at = ext_x , labels = strftime(ext_x,'%H:%M'), tcl = tick_len)
axis(side = 2, at = ext_y, labels = ylab, tcl = tick_len)
axis(side = 3, at = ext_x, labels = NA, tcl = tick_len)
axis(side = 4, at = ext_y, labels = NA, tcl = tick_len)
x_txt <- (xlim[2] - xlim[1])*prc_x+xlim[1]
y_txt <- ylim[2]-(ylim[2] - ylim[1])*prc_y
text(x = x_txt, y_txt,labels = panel.txt)
}
add_night <- function(xlim, ylim){
rise.set = sun.rise.set(xlim[1]+43200, lat) # add mid-day
polygon(x = c(xlim[1], rise.set$sunrise, rise.set$sunrise, xlim[1]), y = c(ylim[1],ylim[1],ylim[2],ylim[2]), col = night_col,border = NA)
polygon(x = c(xlim[2], rise.set$sunset, rise.set$sunset, xlim[2]), y = c(ylim[1],ylim[1],ylim[2],ylim[2]), col = night_col,border = NA)
}
moving_ave <- function(df, window = 18){
out <- df[,2]*NA
l_win <- floor(window/2)
r_win <- window-l_win
end_win <- nrow(df)-r_win
for (i in l_win:end_win){
strt <- i-l_win+1
out[i] <- mean(df[strt:(strt+r_win), 2])
}
df[,2] <- out
return(df)
}
add_models <- function(models){
.empty = sapply(X = models, FUN = function(x) {
lines(moving_ave(x$data), col=x$col, lty = x$lty, lwd = x$lwd)
})
}
add_legend <- function(models, xlim, ylim, prc_x = 0.26, prc_y = 0.06){
y_strt <- ylim[2]-(ylim[2] - ylim[1])*prc_y
y_spc <- (ylim[2] - ylim[1])*0.05
x_len <- (xlim[2] - xlim[1])*0.16
x <- c((xlim[2] - xlim[1])*prc_x+xlim[1], (xlim[2] - xlim[1])*prc_x+xlim[1] + x_len)
for (i in 1:length(models)){
y = y_strt-(i-1)*y_spc
lines(x, c(y,y),
col =models[[i]]$col,
lty = models[[i]]$lty,
lwd = models[[i]]$lwd)
text(x[2],y, models[[i]]$name, pos = 4, cex = 0.65)
}
}
width = 3.37 # single column width for journal
night_col = 'grey90'
height = 1.9
l_mar = 0.35
b_mar = 0.3
t_mar = 0.05
r_mar= 0.05
gapper = 0.15 # space between panels
ylim = c(0,7.55)
xlim = as.POSIXct(c('2009-07-02 22:00', '2009-07-04 02:00', '2009-07-08 22:00', '2009-07-10 02:00'))
default_par = par(no.readonly = TRUE)
#Create plot and save in temporary directory
png(file.path(tempdir(), 'k600_figure.png'), res=300, width=width, height=height, units = 'in')
#
layout(matrix(c(rep(1,10),rep(2,9)),nrow=1)) # 55% on the left panel
par(mai=c(b_mar,l_mar,t_mar,0), omi = c(0,0,0,r_mar),xpd=FALSE,
mgp = c(1.15,.05,0))
plot(c(0,NA),c(0,NA), type='l',
axes = FALSE,
xaxs = 'i', yaxs = 'i',
ylim=ylim,
ylab=expression(k[600]~(m~day^-1)),
xlab=strftime(mean(xlim[1:2]), '%d %b %Y'),
xlim=xlim[1:2])
add_night(xlim[1:2], ylim)
add_models(models)
add_legend(models, xlim, ylim)
add_axes(xlim[1:2], ylim, panel.txt = 'a)')
par(mai=c(b_mar,gapper,t_mar,0))
plot(c(0,NA),c(0,NA), type='l',
axes = FALSE,
xaxs = 'i', yaxs = 'i',
ylim=ylim,
ylab=NA,
xlab=strftime(mean(xlim[3:4]), '%d %b %Y'),
xlim=xlim[3:4])
add_night(xlim[3:4], ylim)
add_models(models)
add_axes(xlim[3:4], ylim, ylabel = NA, panel.txt = 'b)')
#par(xpd=TRUE)
dev.off()
par(default_par)
#quick summary of k600 values for results section
cat('Cole:', mean(k600_cole[,2]), '\n')
cat('Crusius:', mean(k600_crusius[,2]), '\n')
cat('MacIntyre:', mean(k600_macIntyre[,2]), '\n')
cat('Read:', mean(k600_read[,2]), '\n')
cat('Read_soloviev:', mean(k600_soloviev[,2]), '\n')
cat('heiskanen:', mean(k600_heiskanen[,2]), '\n')
cat('vachon:', mean(k600_vachon[,2]), '\n')
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.