# set parameter to prevent overlapping plot
pd = position_dodge(0.1)
## temp. variations
pf_ch = rbind(pf_c,pf_h)
pf_ch$datetime = as.POSIXct(pf_ch$datetime, tz = "Europe/Berlin", format = "%Y-%m-%d %H:%M")
filter1 = pf_ch$datetime < as.POSIXct("2020-10-6 00:00:00")
pf_ch = pf_ch[filter1, ]
filter2 = rect_range$start < as.POSIXct("2020-10-6 00:00:00")
rect_range = rect_range[filter2, ]
p.ch = ggplot(pf_ch, aes(x=datetime, y=Temperature, colour=tnd))+
geom_rect(data = rect_range, mapping=aes(xmin = start,xmax = end,ymin = -Inf, ymax = Inf),
fill = "#3399FF", color = 'NA', alpha = 0.2,inherit.aes = FALSE)+
geom_line(size=1.2)+
xlab("Measuring time") +
ylab(expression(bold(paste('Temperature (',~degree,'C)',sep='')))) +
ylim(15, 40) +
scale_x_datetime(labels = date_format("%m/%d/%y"),
breaks = date_breaks("1 day")) +
scale_color_manual(name = 'Measurements',
breaks = c('control_surface',
'heating_surface'),
values = c("#3399FF","red"),
labels = c('Ambient - sediment surface',
'Heating - sediment surface'))+
theme(
panel.background = element_blank(),
panel.border = element_rect(color = "black", fill=NA, size=2),
legend.position = 'bottom',
legend.title = element_text(size=11, face="bold", color = "black"),
legend.text = element_text(size=10, face="bold"),
axis.title.x = element_text(size = 12, face="bold"),
axis.title.y = element_text(size = 12, face="bold"),
axis.text.x = element_text(size = 10, face="bold"),
axis.text.y = element_text(size = 10, face="bold"))+
guides(color = guide_legend(override.aes = list(size = 1.5)))
p.ch
ggsave('p.ch_temp.png', p.ch, units = 'in', width = 10, height = 4)
lh_temp = pf_ch[pf_ch$Treatment == 'heating', ] %>%
filter(hour(datetime) >= 10 & hour(datetime) <= 14)
View(lh_temp)
mean(lh_temp$Temperature)
min(lh_temp$Temperature)
max(lh_temp$Temperature)
mean(lh_temp$Temperature)
lh_temp = pf_ch[pf_ch$Treatment == 'heating', ] %>%
filter(hour(datetime) >= 11 & hour(datetime) <= 13)
mean(lh_temp$Temperature)
min(lh_temp$Temperature)
max(lh_temp$Temperature)
# packages
library(tidyr)
library(gtools)
library(ggplot2)
library(ggpubr)
library(tidyr)
# Define a function to summarizes data.
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else       length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N    = length2(xx[[col]], na.rm=na.rm),
mean = mean   (xx[[col]], na.rm=na.rm),
sd   = sd     (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
datac <- rename(datac, c("mean" = measurevar))
datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- qt(conf.interval/2 + .5, datac$N-1)
datac$ci <- datac$se * ciMult
return(datac)
}
# define function for maximal change
max_diff = function(data) {
data_max = aggregate(data[, 8:12], by = list(data$day_tank), max)
data_min = aggregate(data[, 8:12], by = list(data$day_tank), min)
colnames(data_max) = c('date_tank', 'cc1_max', 'cc2_max', 'cc3_max', 'cc4_max', 'cc5_max')
colnames(data_min) = c('date_tank', 'cc1_min', 'cc2_min', 'cc3_min', 'cc4_min', 'cc5_min')
data_merge = merge(data_max, data_min)
data_merge$cc1_diff = data_merge$cc1_max - data_merge$cc1_min
data_merge$cc2_diff = data_merge$cc2_max - data_merge$cc2_min
data_merge$cc3_diff = data_merge$cc3_max - data_merge$cc3_min
data_merge$cc4_diff = data_merge$cc4_max -  data_merge$cc4_min
data_merge$cc5_diff = data_merge$cc5_max -  data_merge$cc5_min
date = unique(data$date)
data_merge$date = rep(date, each=2)
data_merge[, c(1,12:17)]
}
depth.hf = read.csv('hf_deltaD.csv', header = T, stringsAsFactors=FALSE)
hf.h5 = subset(depth.hf, phases == 'heating5' & day != 'D28', select = -c(X,line_length))
# packages
library(tidyr)
library(gtools)
library(ggplot2)
library(ggpubr)
library(tidyr)
# Define a function to summarizes data.
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else       length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N    = length2(xx[[col]], na.rm=na.rm),
mean = mean   (xx[[col]], na.rm=na.rm),
sd   = sd     (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
datac <- rename(datac, c("mean" = measurevar))
datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- qt(conf.interval/2 + .5, datac$N-1)
datac$ci <- datac$se * ciMult
return(datac)
}
# define function for maximal change
max_diff = function(data) {
data_max = aggregate(data[, 8:12], by = list(data$day_tank), max)
data_min = aggregate(data[, 8:12], by = list(data$day_tank), min)
colnames(data_max) = c('date_tank', 'cc1_max', 'cc2_max', 'cc3_max', 'cc4_max', 'cc5_max')
colnames(data_min) = c('date_tank', 'cc1_min', 'cc2_min', 'cc3_min', 'cc4_min', 'cc5_min')
data_merge = merge(data_max, data_min)
data_merge$cc1_diff = data_merge$cc1_max - data_merge$cc1_min
data_merge$cc2_diff = data_merge$cc2_max - data_merge$cc2_min
data_merge$cc3_diff = data_merge$cc3_max - data_merge$cc3_min
data_merge$cc4_diff = data_merge$cc4_max -  data_merge$cc4_min
data_merge$cc5_diff = data_merge$cc5_max -  data_merge$cc5_min
date = unique(data$date)
data_merge$date = rep(date, each=2)
data_merge[, c(1,12:17)]
}
depth.hf = read.csv('hf_deltaD.csv', header = T, stringsAsFactors=FALSE)
View(depth.hf)
depth.hf = read.csv('hf_deltaD.csv', header = T, stringsAsFactors=FALSE)
hf.h5 = subset(depth.hf, phases == 'heating5' & day != 'D28', select = -c(X_,line_length))
hf.h5 = spread(data = hf.h5, key = individual, value = deltaL)
depth.lf = read.csv('lf_deltaD.csv', header = T, stringsAsFactors=FALSE)
lf.h3 = subset(depth.lf, phases == 'heating3' & day != 'D28', select = -c(X_,line_length))
View(depth.lf)
depth.hf = read.csv('hf_deltaD.csv', header = T, stringsAsFactors=FALSE)
hf.h5 = subset(depth.hf, phases == 'heating5' & day != 'D28', select = -c(X_,line_length))
hf.h5 = spread(data = hf.h5, key = individual, value = deltaL)
depth.lf = read.csv('lf_deltaD.csv', header = T, stringsAsFactors=FALSE)
lf.h3 = subset(depth.lf, phases == 'heating3' & day != 'D28', select = -c(X,line_length))
lf.h3 = spread(data = lf.h3, key = individual, value = deltaL)
# subset for control and heating
hf.h5.ctrl = subset(hf.h5, treatment == 'control')
hf.h5.ctrl$day_tank = paste0(hf.h5.ctrl$day, '_', hf.h5.ctrl$tank)
hf.h5.heating = subset(hf.h5, treatment == 'heating')
hf.h5.heating$day_tank = paste0(hf.h5.heating$day, '_', hf.h5.heating$tank)
lf.h3.ctrl = subset(lf.h3, treatment == 'control')
lf.h3.ctrl$day_tank = paste0(lf.h3.ctrl$day, '_', lf.h3.ctrl$tank)
lf.h3.heating = subset(lf.h3, treatment == 'heating')
lf.h3.heating$day_tank = paste0(lf.h3.heating$day, '_', lf.h3.heating$tank)
# calculate maximal difference
## for high frequency
h5.c.diff = max_diff(hf.h5.ctrl)
h5.c.diff$treatment = rep('control', 4)
h5.c.diff$frequency = rep('high', 4)
h5.h.diff = max_diff(hf.h5.heating)
h5.h.diff$treatment = rep('heating', 4)
h5.h.diff$frequency = rep('high', 4)
## for low frequency
l3.c.diff = max_diff(lf.h3.ctrl)
l3.c.diff$treatment = rep('control', 4)
l3.c.diff$frequency = rep('low', 4)
l3.h.diff = max_diff(lf.h3.heating)
l3.h.diff$treatment = rep('heating', 4)
l3.h.diff$frequency = rep('low', 4)
## merge data
al.diff.raw = rbind(h5.c.diff, h5.h.diff, l3.h.diff)
al.diff = gather(al.diff.raw, key = individual, value = maximal_change, 2:6, na.rm = T)
write.csv(al.diff, 'al.diff.csv', row.names = F)
# loading data
lumino = read.csv('lumino.csv', header =T)
lumino.lhf = subset(lumino, Day == 'H3' | Day == 'L3' | Day == 'M1'| Day == 'M2'| Day == 'M3')
lumino.lhf = subset(lumino.lhf, LnB == 'B')
al.diff = read.csv('al.diff.csv', header =T)
# subset and organise
lumino.c = subset(lumino.lhf, Setting == 'monitoring')
lumino.c = subset(lumino.c, Treatment != 'heating')
lumino.l = subset(lumino.lhf, Setting == 'low')
lumino.h = subset(lumino.lhf, Setting == 'high')
lumino.c$tank_rep = paste0(lumino.c$Tank, '_', lumino.c$Beplicat)
lumino.l$tank_rep = paste0(lumino.l$Tank, '_', lumino.l$Beplicat)
lumino.h$tank_rep = paste0(lumino.h$Tank, '_', lumino.h$Beplicat)
lumino.clh = rbind(lumino.c, lumino.h, lumino.l)
# give tank and individule name
lumino.clh$Tank[lumino.clh$Tank == 'C1'] = 'tank1'
lumino.clh$Tank[lumino.clh$Tank == 'C2'] = 'tank2'
lumino.clh$Tank[lumino.clh$Tank == 'H2'] = 'tank3'
lumino.clh$Tank[lumino.clh$Tank == 'H1'] = 'tank4'
al.diff = separate(al.diff, col = 'date_tank', into = c('date', 'tank'))
al.diff$individual[al.diff$individual == 'cc1_diff'] = 'B1'
al.diff$individual[al.diff$individual == 'cc2_diff'] = 'B2'
al.diff$individual[al.diff$individual == 'cc3_diff'] = 'B3'
al.diff$individual[al.diff$individual == 'cc4_diff'] = 'B4'
al.diff$individual[al.diff$individual == 'cc5_diff'] = 'B5'
#summarize data for lumino and depth
## creat labels
al.diff$frequency_treatment = paste0(al.diff$frequency, '_', al.diff$treatment)
lumino.clh$setting_tank_individual_treatment = paste0(lumino.clh$Setting, '_', lumino.clh$Tank, '_', lumino.clh$Beplicat, '_', lumino.clh$Treatment)
lumino.clh = subset(lumino.clh, area_areahole_blank >= 0 & Depth != 0.5 & Depth != 4.0 & Depth != 5.0)
rownames(lumino.clh) = 1:100
## summarize data
diff_mean = aggregate(al.diff$maximal_change, by = list(Category = al.diff$frequency_treatment), FUN = mean)
colnames(diff_mean)[2] = 'maximal_change'
diff_mean[,1] = c('control', 'high', 'low')
lumino_sum = aggregate(lumino.clh$area_adjusted_perc, list(Category = lumino.clh$setting_tank_individual_treatment), FUN=sum)
colnames(lumino_sum)[2] = 'bioturbated_area'
lumino_sum$bioturbated_area = (lumino_sum$bioturbated_area)/100 * 35.23865 *5
lumino_sum = separate(lumino_sum, col = 'Category', into = c('setting', 'tank', 'rep', 'treatment'))
lumino_mean = aggregate(lumino_sum$bioturbated_area, list(Category = lumino_sum$setting), FUN=mean)
colnames(lumino_mean)[2] = 'bioturbated_area'
lumino_mean[,1] = c('high', 'low', 'control')
lumino_diff_mean = merge(diff_mean, lumino_mean)
colnames(lumino_diff_mean)[1] = 'setting'
lumino_diff_mean$survival_rate= c(100, 80, 45)
lumino_diff_mean$setting = c('Control', 'High frequency', 'Low frequency')
ldm = lumino_diff_mean
ldm$mortality_rate = 100 - ldm$survival_rate
ldm.points = ldm[, c(1,2,5)]
p.ldm = ggplot(data = ldm, aes(y = maximal_change, x = mortality_rate,
color = setting, size = setting)) +
geom_point() +
geom_point(data = ldm.points, aes(y = maximal_change, x = mortality_rate),
color = 'black', size = 3.5, shape = 10)+
coord_cartesian(xlim = c(-10, 75), ylim = c(-0.2, 5)) +
scale_size_manual(name = 'Bioturbated area',
breaks = c("Control", "High frequency", "Low frequency"),
values = ldm$bioturbated_area/1.8, guide = F) +
scale_colour_manual(name = 'Treatments',
values = c('#3399FF', '#FF9900', '#FF0000'),
breaks = c("Control", "High frequency", "Low frequency"),
labels = c("Control ambient temperature",
"3-day cycle heatwaves",
"6-day cycle heatwaves"))+
ylab('average maximal position change (mm)')+
xlab('mortality rate (%)') +
theme(
plot.title = element_blank(),
panel.background = element_blank(),
legend.title = element_text(size=24, face="bold", color = "black"),
legend.text = element_text(size=24, face="bold"),
legend.position=c(0.3, 0.14),
legend.background = element_rect(fill=alpha('white', 0.8),
size=0.5, linetype="solid",
colour ="NA"),
legend.key = element_blank(),
axis.title.x = element_text(size = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold"),
axis.text.y = element_text(size = 28, face="bold"),
axis.ticks = element_line(size = 2),
axis.ticks.length = unit(0.08, "inch"),
plot.margin = unit(c(0.1,0.1,0.1,0.1), "in"))+
guides(color = guide_legend(override.aes = list(size = 5)))
p.ldm = p.ldm +
geom_segment(aes(x=0, xend=60, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=5, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
p.ldm
ggsave('ldm.png', p.ldm, units = 'in', width = 12, height = 9,dpi=600)
# p.ldm = ggplot(data = ldm, aes(x = maximal_change, y = survival_rate, color = setting)) +
#   geom_point(aes(size = setting)) +
#   scale_size_manual(name = 'Bioturbated area',
#                     breaks = c("Control", "High frequency", "Low frequency"),
#                     values = ldm$bioturbated_area*2, guide = F) +
#   scale_colour_manual(name = 'Heat wave frequency',
#                       values = c('#3399FF', '#FF9900', '#FF0000'),
#                       breaks = c("Control", "High frequency", "Low frequency"),
#                       labels = c("Control", "High frequency", "Low frequency"))+
#   xlab('average maximal position change (mm)')+
#   ylab('survival rate (%)') +
#   xlim(1.75, 4.25) +
#   ylim(40, 110) +
#   theme_bw()+
#   theme(legend.position='bottom',
#         legend.title = element_text(size=12, face="bold", color = "blue"),
#         legend.text = element_text(size=11, face="bold"),
#         strip.text.x = element_text(size = 14, face = "bold"),
#         axis.title.x = element_text(size = 15, face="bold"),
#         axis.title.y = element_text(size = 15, face="bold"),
#         axis.text.x = element_text(size = 12, face="bold"),
#         axis.text.y = element_text(size = 12, face="bold"))+
#   guides(colour = guide_legend(override.aes = list(size=4)))
# p.ldm
# ggsave('ldm.png', p.ldm, units = 'in', width = 7.5, height = 7,dpi=600)
View(ldm)
p.ldm = ggplot(data = ldm, aes(y = bioturbated_area, x = setting,
color = setting)) +
geom_point()
p.ldm
p.ldm = ggplot(data = ldm, aes(y = bioturbated_area/5, x = setting,
color = setting)) +
geom_point()+
scale_colour_manual(name = 'Treatments',
values = c('#3399FF', '#FF9900', '#FF0000'),
breaks = c("Control", "High frequency", "Low frequency"),
labels = c("Control ambient temperature",
"3-day cycle heatwaves",
"6-day cycle heatwaves"))+
ylab('total bioturbbation area')+
xlab('heatwave scenarios (%)') +
theme(
plot.title = element_blank(),
panel.background = element_blank(),
legend.title = element_text(size=24, face="bold", color = "black"),
legend.text = element_text(size=24, face="bold"),
legend.position=c(0.3, 0.14),
legend.background = element_rect(fill=alpha('white', 0.8),
size=0.5, linetype="solid",
colour ="NA"),
legend.key = element_blank(),
axis.title.x = element_text(size = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold"),
axis.text.y = element_text(size = 28, face="bold"),
axis.ticks = element_line(size = 2),
axis.ticks.length = unit(0.08, "inch"),
plot.margin = unit(c(0.1,0.1,0.1,0.1), "in"))+
guides(color = guide_legend(override.aes = list(size = 5)))
p.ldm
p.ldm = p.ldm +
geom_segment(aes(x=0, xend=60, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=5, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
p.ldm
p.ldm = p.ldm +
geom_segment(aes(x=0, xend=3, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=30, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
p.ldm
p.ldm = p.ldm +
geom_segment(aes(x=0, xend=1, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=30, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
p.ldm
p.ldm = ggplot(data = ldm, aes(y = bioturbated_area/5, x = setting,
color = setting)) +
geom_point()+
scale_colour_manual(name = 'Treatments',
values = c('#3399FF', '#FF9900', '#FF0000'),
breaks = c("Control", "High frequency", "Low frequency"),
labels = c("Control ambient temperature",
"3-day cycle heatwaves",
"6-day cycle heatwaves"))+
ylab('total bioturbbation area')+
xlab('heatwave scenarios (%)') +
theme(
plot.title = element_blank(),
panel.background = element_blank(),
legend.title = element_text(size=24, face="bold", color = "black"),
legend.text = element_text(size=24, face="bold"),
legend.position=c(0.3, 0.14),
legend.background = element_rect(fill=alpha('white', 0.8),
size=0.5, linetype="solid",
colour ="NA"),
legend.key = element_blank(),
axis.title.x = element_text(size = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold"),
axis.text.y = element_text(size = 28, face="bold"),
axis.ticks = element_line(size = 2),
axis.ticks.length = unit(0.08, "inch"),
plot.margin = unit(c(0.1,0.1,0.1,0.1), "in"))+
guides(color = guide_legend(override.aes = list(size = 5)))
p.ldm = p.ldm +
geom_segment(aes(x=0, xend=1, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=30, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
p.ldm
p.ldm = ggplot(data = ldm, aes(y = bioturbated_area/5, x = setting,
color = setting)) +
geom_point(size = 3.5, shape = 10)+
scale_colour_manual(name = 'Treatments',
values = c('#3399FF', '#FF9900', '#FF0000'),
breaks = c("Control", "High frequency", "Low frequency"),
labels = c("Control ambient temperature",
"3-day cycle heatwaves",
"6-day cycle heatwaves"))+
ylab('total bioturbbation area')+
xlab('heatwave scenarios (%)') +
theme(
plot.title = element_blank(),
panel.background = element_blank(),
legend.title = element_text(size=24, face="bold", color = "black"),
legend.text = element_text(size=24, face="bold"),
legend.position=c(0.3, 0.14),
legend.background = element_rect(fill=alpha('white', 0.8),
size=0.5, linetype="solid",
colour ="NA"),
legend.key = element_blank(),
axis.title.x = element_text(size = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold"),
axis.text.y = element_text(size = 28, face="bold"),
axis.ticks = element_line(size = 2),
axis.ticks.length = unit(0.08, "inch"),
plot.margin = unit(c(0.1,0.1,0.1,0.1), "in"))+
guides(color = guide_legend(override.aes = list(size = 5)))
p.ldm = p.ldm +
geom_segment(aes(x=0, xend=3, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=30, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
p.ldm
################# Plot temperature data #################
######loading packages######
library(dplyr)
library(ggpubr)
library(ggplot2)
library(gridExtra)
library(scales)
library(cowplot)
library(tidyr)
library(lubridate)
######loading data######
pf = read.csv('sorted_pf1.csv', header = T, row.names = 1)
pf$time = substr(as.POSIXct(sprintf("%04.0f", pf$time), format='%H%M'), 12, 16)
pf$minute  = factor(pf$time)
pf$minute = strptime(pf$minute, "%H:%M")
pf$minute = format(pf$minute, "%M")
pf$minute = as.numeric(pf$minute)
pf = subset(pf, Depth != '6cm'& Depth != '3cm')
#pf = subset(pf, Tank == 'tank2'| Tank == 'tank4')
pf$tnd = paste(pf$Treatment, "_", pf$Depth,sep='')
pf = subset(pf, minute %% 5 == 0)
rect_range = read.csv('tide_schedule.csv', header = T)
rect_range$start = as.POSIXct(rect_range$start, tz = "Europe/Berlin", format = "%m/%d/%Y %H:%M")
rect_range$end = as.POSIXct(rect_range$end, tz = "Europe/Berlin", format = "%m/%d/%Y %H:%M")
#treatments
pf_control = subset(pf, select = -Sensor, Treatment == 'control')
pf_heat = subset(pf, select = -Sensor, Treatment == 'heating')
pf_c = pf_control %>%
group_by(datetime, Treatment, Depth,tnd) %>%
summarise_all("mean")
pf_h = pf_heat %>%
group_by(datetime, Treatment, Depth,tnd) %>%
summarise_all("mean")
######plotting######
# set parameter to prevent overlapping plot
pd = position_dodge(0.1)
## temp. variations
pf_ch = rbind(pf_c,pf_h)
pf_ch$datetime = as.POSIXct(pf_ch$datetime, tz = "Europe/Berlin", format = "%Y-%m-%d %H:%M")
filter1 = pf_ch$datetime < as.POSIXct("2020-10-6 00:00:00")
pf_ch = pf_ch[filter1, ]
filter2 = rect_range$start < as.POSIXct("2020-10-6 00:00:00")
rect_range = rect_range[filter2, ]
p.ch = ggplot(pf_ch, aes(x=datetime, y=Temperature, colour=tnd))+
geom_rect(data = rect_range, mapping=aes(xmin = start,xmax = end,ymin = -Inf, ymax = Inf),
fill = "#3399FF", color = 'NA', alpha = 0.2,inherit.aes = FALSE)+
geom_line(size=1.2)+
xlab("Measuring time") +
ylab(expression(bold(paste('Temperature (',~degree,'C)',sep='')))) +
ylim(15, 40) +
scale_x_datetime(labels = date_format("%m/%d/%y"),
breaks = date_breaks("1 day")) +
scale_color_manual(name = 'Measurements',
breaks = c('control_surface',
'heating_surface'),
values = c("#3399FF","red"),
labels = c('Ambient - sediment surface',
'Heating - sediment surface'))+
theme(
panel.background = element_blank(),
panel.border = element_rect(color = "black", fill=NA, size=2),
legend.position = 'bottom',
legend.title = element_text(size=11, face="bold", color = "black"),
legend.text = element_text(size=10, face="bold"),
axis.title.x = element_text(size = 12, face="bold"),
axis.title.y = element_text(size = 12, face="bold"),
axis.text.x = element_text(size = 10, face="bold"),
axis.text.y = element_text(size = 10, face="bold"))+
guides(color = guide_legend(override.aes = list(size = 1.5)))
p.ch
ggsave('p.ch_temp.png', p.ch, units = 'in', width = 10, height = 4)
lh_temp = pf_ch[pf_ch$Treatment == 'heating', ] %>%
filter(hour(datetime) >= 11 & hour(datetime) <= 13)
min(lh_temp$Temperature) # 28.1825
max(lh_temp$Temperature) # 37.4725
mean(lh_temp$Temperature)
max(lh_temp$Temperature)
min(lh_temp$Temperature)
(28+37)/2
View(pf_ch)
lc_temp = pf_ch[pf_ch$Treatment == 'control', ] %>%
filter(hour(datetime) >= 11 & hour(datetime) <= 13)
min(lh_temp$Temperature) # 28.1825
min(lc_temp$Temperature) # 28.1825
max(lc_temp$Temperature) # 37.4725
mean(lc_temp$Temperature) #32.22539
sd(lh_temp$Temperature)
sd(lc_temp$Temperature) #
