# ------------------------------------------------------------- # # Compare the datasets in Klepstand monitor project # Voordelta (winter and summer) and betonbak # # # ------------------------------------------------------------- #load necessary packages require(forecast) require(lmtest) #grangertest() function require(stringr) #load the files Perc_zo = read.csv("Voordelta_zo_Perc.csv") Perc_wi = read.csv("Voordelta_wi_Perc.csv") Flow_wi = read.csv("Flow_wi.csv") Flow_zo = read.csv("Flow_zo.csv") Chla_wi = read.csv("chla_wi.csv") Chla_zo = read.csv("chla_zo.csv") # Betonbak datasets Perc_BB = read.csv("Betonbak_Perc.csv") Chla_BB = read.csv("chla_betonbak.csv") # --------------------------------------------------------------------------------------------------- # RWS tidal data Beton_tide = read.csv2("./Raw datasets/NVT_WATHTBRKD_YERSKE.csv", skip = 5, head = FALSE)[ ,-6] # last column is empty # Give columns more meaningful names colnames(Beton_tide) = c("Date", "Time", "Details", "Location", "Tide") Beton_tide$units = str_sub(Beton_tide$Details,-2,-1) plot(Beton_tide$Tide) Beton_tide$Date = as.Date(Beton_tide$Date, format = "%d/%m/%Y") Beton_tide$Year = str_sub(Beton_tide$Date, 1, 4) Beton_tide$Month = str_sub(Beton_tide$Date, 6, 7) Beton_tide$Day = str_sub(Beton_tide$Date, 9, 10) Beton_tide$Hour = str_sub(Beton_tide$Time, 1, 2) Beton_tide$Minute = str_sub(Beton_tide$Time, 4, 5) Beton_tide$Second = str_sub(Beton_tide$Time, 7, 8) attach(Beton_tide) Beton_tide$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Beton_tide$Datetime = as.POSIXct(Beton_tide$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Beton_tide$Datetime) # the ordering of the df is reversed. fix by: Beton_tide = Beton_tide[order(nrow(Beton_tide):1),] # Match the time range with the betonbak range of 17 days Beton_tide = Beton_tide[Beton_tide$Datetime > "2020-08-12 11:59:50" & Beton_tide$Datetime < "2020-08-29 11:59:50", ] with(Beton_tide, difftime(max(Beton_tide$Datetime, na.rm = TRUE), min(Beton_tide$Datetime, na.rm = TRUE))) # ------------------------------------------------------------------------------------------------------------ # standardize all datasets to 30-minute intervals (not moving average of 30 minutes, which retains all points) # ------------------------------------------------------------------------------------------------------------ #dates in flow and perc files are converting incorrectly; remake the date column attach(Flow_wi) Flow_wi$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Flow_wi$Datetime = as.POSIXct(Flow_wi$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") attach(Flow_zo) Flow_zo$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Flow_zo$Datetime = as.POSIXct(Flow_zo$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Flow_zo$Datetime) attach(Perc_zo) Perc_zo$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Perc_zo$Datetime = as.POSIXct(Perc_zo$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Perc_zo$Datetime) attach(Perc_wi) Perc_wi$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Perc_wi$Datetime = as.POSIXct(Perc_wi$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Perc_wi$Datetime) attach(Perc_BB) Perc_BB$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Perc_BB$Datetime = as.POSIXct(Perc_BB$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Perc_BB$Datetime) # --------------------------------------------------------------------------------------------------- # Tide data for Voordelta (Roompot Buiten buoy) for spring / autumn Autumn_tide = read.csv2("./Raw datasets/NVT_WATHTBRKD_ROOMPBTN_autumn.csv", skip = 5, head = FALSE)[ ,-6] # Give columns more meaningful names colnames(Autumn_tide) = c("Date", "Time", "Details", "Location", "Tide") Autumn_tide$units = str_sub(Autumn_tide$Details,-2,-1) plot(Autumn_tide$Tide) Autumn_tide$Date = as.Date(Autumn_tide$Date, format = "%d/%m/%Y") Autumn_tide$Year = str_sub(Autumn_tide$Date, 1, 4) Autumn_tide$Month = str_sub(Autumn_tide$Date, 6, 7) Autumn_tide$Day = str_sub(Autumn_tide$Date, 9, 10) Autumn_tide$Hour = str_sub(Autumn_tide$Time, 1, 2) Autumn_tide$Minute = str_sub(Autumn_tide$Time, 4, 5) Autumn_tide$Second = str_sub(Autumn_tide$Time, 7, 8) attach(Autumn_tide) Autumn_tide$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Autumn_tide$Datetime = as.POSIXct(Autumn_tide$Datetime,format='%Y/%m/%d %H:%M:%S',tz="CET") class(Autumn_tide$Datetime) # the ordering of the df is reversed. fix by: Autumn_tide = Autumn_tide[order(nrow(Autumn_tide):1),] # Match the time range with the autumn voordelta range of 17 days Autumn_tide = Autumn_tide[Autumn_tide$Datetime > min(Perc_wi$Datetime) & Autumn_tide$Datetime < max(Perc_wi$Datetime), ] with(Autumn_tide, difftime(max(Autumn_tide$Datetime, na.rm = TRUE), min(Autumn_tide$Datetime, na.rm = TRUE))) Autumn_tide$Tide = as.numeric(as.character(Autumn_tide$Tide)) Autumn_tide$Tide = Autumn_tide$Tide/100 par(mar = c(6,6,2,6)) plot(ma(Flow_wi$Speed, order = 3, centre = TRUE) , type = "l", las = 1, ylab = "", ylim = c(0, 1), xaxt = "n", ann = FALSE) mtext(side = 2, line = 3, expression("Flow velocity (m s"^{-1}*")")) mtext(side = 3, "Autumn - Voordelta") par(new=TRUE) plot(Autumn_tide$Tide, ylim = c(-4, 4), type = "l", col = "red", axes = F, xlab = NA, ylab = NA, lty = 3) mtext(side = 4, line = 3, "Tidal level (m)") axis(side = 4, at = seq(from = min(Autumn_tide$Tide), to = max(Autumn_tide$Tide), length.out = 5), labels = round(seq(from = min(Autumn_tide$Tide), to = max(Autumn_tide$Tide), length.out = 5), 0), las = 1) legend("topleft", c("Flow", "Tide"),bty="n",cex=1.25, lty = c(1, 3), col = c("black", "red"), lwd = c(2,2)) axis(side = 1, at = seq(from = 0, to = 2448, by = 144), labels = c(0:17)) mtext(side = 1, line = 3.5, "Time (Days)") Min_30 = cut(Autumn_tide$Datetime, breaks = "30 min") Autumn_tide_30 = aggregate(Tide ~ Min_30, Autumn_tide, mean) # --------------------------------------------------------------------------------------------------- Spring_tide = read.csv2("./Raw datasets/NVT_WATHTBRKD_ROOMPBTN_spring.csv", skip = 5, head = FALSE)[ ,-6] # Give columns more meaningful names colnames(Spring_tide) = c("Date", "Time", "Details", "Location", "Tide") Spring_tide$units = str_sub(Spring_tide$Details,-2,-1) plot(Spring_tide$Tide) Spring_tide$Date = as.Date(Spring_tide$Date, format = "%d/%m/%Y") Spring_tide$Year = str_sub(Spring_tide$Date, 1, 4) Spring_tide$Month = str_sub(Spring_tide$Date, 6, 7) Spring_tide$Day = str_sub(Spring_tide$Date, 9, 10) Spring_tide$Hour = str_sub(Spring_tide$Time, 1, 2) Spring_tide$Minute = str_sub(Spring_tide$Time, 4, 5) Spring_tide$Second = str_sub(Spring_tide$Time, 7, 8) attach(Spring_tide) Spring_tide$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Spring_tide$Datetime = as.POSIXct(Spring_tide$Datetime,format='%Y/%m/%d %H:%M:%S',tz="CET") class(Spring_tide$Datetime) # the ordering of the df is reversed. fix by: Spring_tide = Spring_tide[order(nrow(Spring_tide):1),] # Match the time range with the autumn voordelta range of 17 days Spring_tide = Spring_tide[Spring_tide$Datetime > min(Perc_zo$Datetime) & Spring_tide$Datetime < max(Perc_zo$Datetime), ] with(Spring_tide, difftime(max(Spring_tide$Datetime, na.rm = TRUE), min(Spring_tide$Datetime, na.rm = TRUE))) Spring_tide$Tide = as.numeric(as.character(Spring_tide$Tide)) Spring_tide$Tide = Spring_tide$Tide/100 # 30-minute averaged: tide vs. flow (regression curves) par(mar = c(6,6,2,6)) plot(ma(Flow_zo$Speed, order = 3, centre = TRUE) , type = "l", las = 1, ylab = "", ylim = c(0, 1), xaxt = "n", ann = FALSE) mtext(side = 2, line = 3, expression("Flow velocity (m s"^{-1}*")")) mtext(side = 3, "Spring - Voordelta") par(new=TRUE) plot(Spring_tide$Tide, ylim = c(-4, 4), type = "l", col = "red", axes = F, xlab = NA, ylab = NA, lty = 3) mtext(side = 4, line = 3, "Tidal level (m)") axis(side = 4, at = seq(from = min(Autumn_tide$Tide), to = max(Autumn_tide$Tide), length.out = 5), labels = round(seq(from = min(Autumn_tide$Tide), to = max(Autumn_tide$Tide), length.out = 5), 0), las = 1) legend("topleft", c("Flow", "Tide"),bty="n",cex=1.25, lty = c(1, 3), col = c("black", "red"), lwd = c(2,2)) axis(side = 1, at = seq(from = 0, to = 2448, by = 144), labels = c(0:17)) mtext(side = 1, line = 3.5, "Time (Days)") Min_30 = cut(Spring_tide$Datetime, breaks = "30 min") Spring_tide_30 = aggregate(Tide ~ Min_30, Spring_tide, mean) # 30-minute averaged: tide vs. flow (regression curves) plot(Flow_zo_30$Speed ~ Spring_tide_30$Tide, ann = FALSE, las = 1) mtext(side = 3, "Spring - Voordelta") mtext(side = 2, line = 3, expression("Flow velocity (m s"^{-1}*")")) mtext(side = 1, line = 3, "Tidal level (m)") model1 <-lm(Flow_zo_30$Speed ~ Spring_tide_30$Tide+I(Spring_tide_30$Tide^2)) plotdata <- cbind(Spring_tide_30$Tide, predict(model1)) lines(plotdata[order(Spring_tide_30$Tide),], col = "blue", lwd = 1.5) mtext(paste("R2 =",round(summary(model1)$r.squared,3)), side=3, line=-1.5, adj=0.05, cex = 1) mtext(paste("beta =",round(model1$coefficients[2],3)), side=3, line=-2.5, adj=0.05, cex = 1) mtext(paste("p =",round(summary(model1)$coef[2,4],5)), side=3, line=-3.5, adj=0.05, cex = 1) # --------------------------------------------------------------------------------------------------- # length of the datasets (should be longer/shorter by #times of the measurement interval) length(ma(Flow_zo$Speed, order = 3, centre = TRUE) ) length(ma(Chla_zo$Chl.a.ug.l., order = 30, centre = TRUE) ) length(ma(Perc_zo$Perc.1, order = 30, centre = TRUE) ) length(ma(Flow_wi$Speed, order = 3, centre = TRUE) ) length(ma(Chla_wi$Chl.a.ug.l., order = 30, centre = TRUE)) length(ma(Perc_wi$Perc.8, order = 30, centre = TRUE) ) length(ma(Chla_BB$Chl.a.ug.l., order = 30, centre = TRUE)) length(ma(Perc_BB$Perc.8, order = 30, centre = TRUE) ) library(anytime) grangertest(ma(Flow_zo$Speed, order = 3, centre = TRUE) , ma(Chla_zo$Chl.a.ug.l., order = 30, centre = TRUE) ) plot(ma(Flow_zo$Speed, order = 3, centre = TRUE) , ma(Chla_zo$Chl.a.ug.l., order = 30, centre = TRUE)) #chl datasets do not have the individual date/time columns: first create them Chla_wi$Year = str_sub(Chla_wi$Date, 1, 4) Chla_wi$Month = str_sub(Chla_wi$Date, 6, 7) Chla_wi$Day = str_sub(Chla_wi$Date, 9, 10) Chla_wi$Hour = str_sub(Chla_wi$Date, 12, 13) Chla_wi$Minute = str_sub(Chla_wi$Date, 15, 16) Chla_wi$Second = str_sub(Chla_wi$Date, 18, 19) attach(Chla_wi) Chla_wi$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Chla_wi$Datetime = as.POSIXct(Chla_wi$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Chla_wi$Datetime) Chla_zo$Year = str_sub(Chla_zo$Date, 1, 4) Chla_zo$Month = str_sub(Chla_zo$Date, 6, 7) Chla_zo$Day = str_sub(Chla_zo$Date, 9, 10) Chla_zo$Hour = str_sub(Chla_zo$Date, 12, 13) Chla_zo$Minute = str_sub(Chla_zo$Date, 15, 16) Chla_zo$Second = str_sub(Chla_zo$Date, 18, 19) attach(Chla_zo) Chla_zo$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Chla_zo$Datetime = as.POSIXct(Chla_zo$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Chla_zo$Datetime) Chla_BB$Year = str_sub(Chla_BB$Date, 1, 4) Chla_BB$Month = str_sub(Chla_BB$Date, 6, 7) Chla_BB$Day = str_sub(Chla_BB$Date, 9, 10) Chla_BB$Hour = str_sub(Chla_BB$Date, 12, 13) Chla_BB$Minute = str_sub(Chla_BB$Date, 15, 16) Chla_BB$Second = str_sub(Chla_BB$Date, 18, 19) attach(Chla_BB) Chla_BB$Datetime = ISOdate((Year),Month,Day,Hour,Minute,Second,tz="CET") Chla_BB$Datetime = as.POSIXct(Chla_BB$Datetime,format='%Y/%M/%d %H:%M:%S',tz="CET") class(Chla_BB$Datetime) #check the date/time ranges are equal between the datasets range(Perc_zo$Datetime) range(Flow_zo$Datetime) range(Chla_zo$Datetime) range(Perc_wi$Datetime) range(Flow_wi$Datetime) #one day too long range(Chla_wi$Datetime) range(Perc_BB$Datetime) range(Chla_BB$Datetime) range(Beton_tide$Datetime) #remove extra day Flow_wi = Flow_wi[Flow_wi$Datetime >= min(Chla_wi$Datetime) & Flow_wi$Datetime <= max(Chla_wi$Datetime), ] # --------------------------------------- # Note that the Perc_zo/wi/BB datasets are 10 seconds ahead # add 10 seconds to match (remove the last datapoint, which is 00:09:50 after the last point for chla/flow) Perc_zo$Datetime = Perc_zo$Datetime + 10 Perc_zo = Perc_zo[Perc_zo$Datetime >= min(Chla_zo$Datetime) & Perc_zo$Datetime <= max(Chla_zo$Datetime), ] Perc_wi$Datetime = Perc_wi$Datetime + 10 Perc_wi = Perc_wi[Perc_wi$Datetime >= min(Chla_wi$Datetime) & Perc_wi$Datetime <= max(Chla_wi$Datetime), ] Perc_BB$Datetime = Perc_BB$Datetime + 10 Perc_BB = Perc_BB[Perc_BB$Datetime >= min(Chla_BB$Datetime) & Perc_BB$Datetime <= max(Chla_BB$Datetime), ] #every set should match now range(Perc_zo$Datetime) range(Flow_zo$Datetime) range(Chla_zo$Datetime) range(Perc_wi$Datetime) range(Flow_wi$Datetime) range(Chla_wi$Datetime) range(Perc_BB$Datetime) range(Chla_BB$Datetime) range(Beton_tide$Datetime) ############################################################################################## # --------------------------------------------------------------------------------- # Average by day for Gape, Temp, Chl-a and turbidity # --------------------------------------------------------------------------------- # ------------------ # summer # ------------------ One_day = cut(Chla_zo$Datetime, breaks = "1440 min") chla_zo_day = aggregate(Chl.a.ug.l. ~ One_day, Chla_zo, mean) Temp_zo_day = aggregate(Temp..deg.C. ~ One_day, Chla_zo, mean) Turb_zo_day = aggregate(Turb...M.FTU. ~ One_day, Chla_zo, mean) # one-day averaging is not meaningful for flow # One_day = cut(Flow_zo$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "1440 min") # Flow_zo_day = aggregate(Speed ~ One_day, Flow_zo, mean) One_day = cut(Perc_zo$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "1440 min") Perc1_zo_day = aggregate(Perc.1 ~ One_day, Perc_zo, mean) Perc2_zo_day = aggregate(Perc.2 ~ One_day, Perc_zo, mean) Perc5_zo_day = aggregate(Perc.5 ~ One_day, Perc_zo, mean) #Perc6_zo = aggregate(Ch6.Perc.6 ~ One_day, Perc_zo, mean) # Perc7_zo_day = aggregate(Perc.7 ~ One_day, Perc_zo, mean) # Perc8_zo_day = aggregate(Perc.8 ~ One_day, Perc_zo, mean) Perc_zo_day = cbind(Perc1_zo_day, Perc2_zo_day[2], Perc5_zo_day[2]) ENVI_zo_day = cbind(chla_zo_day, Temp_zo_day[2], Turb_zo_day[2]) write.csv(Perc_zo_day, "Perc_zo_day.csv") write.csv(ENVI_zo_day, "ENVI_zo_day.csv") # ------------------ # winter # ------------------ One_day = cut(Chla_wi$Datetime, breaks = "1440 min") chla_wi_day = aggregate(Chl.a.ug.l. ~ One_day, Chla_wi, mean) Temp_wi_day = aggregate(Temp..deg.C. ~ One_day, Chla_wi, mean) Turb_wi_day = aggregate(Turb...M.FTU. ~ One_day, Chla_wi, mean) One_day = cut(Perc_wi$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "1440 min") Perc1_wi_day = aggregate(Perc.1 ~ One_day, Perc_wi, mean) Perc6_wi_day = aggregate(Perc.6 ~ One_day, Perc_wi, mean) Perc7_wi_day = aggregate(Perc.7 ~ One_day, Perc_wi, mean) Perc8_wi_day = aggregate(Perc.8 ~ One_day, Perc_wi, mean) Perc_wi_day = cbind(Perc1_wi_day, Perc6_wi_day[2], Perc7_wi_day[2], Perc8_wi_day[2]) ENVI_wi_day = cbind(chla_wi_day, Temp_wi_day[2], Turb_wi_day[2]) write.csv(Perc_wi_day, "Perc_wi_day.csv") write.csv(ENVI_wi_day, "ENVI_wi_day.csv") # ------------------ # Betonbak # ------------------ One_day = cut(Chla_BB$Datetime, breaks = "1440 min") chla_BB_day = aggregate(Chl.a.ug.l. ~ One_day, Chla_BB, mean) Temp_BB_day = aggregate(Temp..deg.C. ~ One_day, Chla_BB, mean) Turb_BB_day = aggregate(Turb...M.FTU. ~ One_day, Chla_BB, mean) One_day = cut(Perc_BB$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "1440 min") Perc1_BB_day = aggregate(Perc.1 ~ One_day, Perc_BB, mean) Perc2_BB_day = aggregate(Perc.2 ~ One_day, Perc_BB, mean) Perc3_BB_day = aggregate(Perc.3 ~ One_day, Perc_BB, mean) Perc4_BB_day = aggregate(Perc.4 ~ One_day, Perc_BB, mean) Perc5_BB_day = aggregate(Perc.5 ~ One_day, Perc_BB, mean) Perc6_BB_day = aggregate(Perc.6 ~ One_day, Perc_BB, mean) Perc7_BB_day = aggregate(Perc.7 ~ One_day, Perc_BB, mean) Perc8_BB_day = aggregate(Perc.8 ~ One_day, Perc_BB, mean) Perc_BB_day = cbind(Perc1_BB_day, Perc2_BB_day[2], Perc3_BB_day[2], Perc4_BB_day[2], Perc5_BB_day[2], Perc6_BB_day[2], Perc7_BB_day[2], Perc8_BB_day[2]) ENVI_BB_day = cbind(chla_BB_day, Temp_BB_day[2], Turb_BB_day[2]) write.csv(Perc_BB_day, "Perc_BB_day.csv") write.csv(ENVI_BB_day, "ENVI_BB_day.csv") ###################################################################################### # --------------------------------------------------------------------------------- # 30 - minute averaged (not rolling average) # Very similar to 30-minute rolling average, only the number of datapoints is # drastically reduced # --------------------------------------------------------------------------------- # summer Min_30 = cut(Chla_zo$Datetime, breaks = "30 min") chla_zo_30 = aggregate(Chl.a.ug.l. ~ Min_30, Chla_zo, mean) Temp_zo_30 = aggregate(Temp..deg.C. ~ Min_30, Chla_zo, mean) Turb_zo_30 = aggregate(Turb...M.FTU. ~ Min_30, Chla_zo, mean) Min_30 = cut(Flow_zo$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "30 min") Flow_zo_30 = aggregate(Speed ~ Min_30, Flow_zo, mean) Min_30 = cut(Perc_zo$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "30 min") Perc1_zo = aggregate(Perc.1 ~ Min_30, Perc_zo, mean) Perc2_zo = aggregate(Perc.2 ~ Min_30, Perc_zo, mean) Perc5_zo = aggregate(Perc.5 ~ Min_30, Perc_zo, mean) Perc_zo_30 = cbind(Perc1_zo, Perc2_zo[2], Perc5_zo[2]) ENVI_zo_30 = cbind(chla_zo_30, Temp_zo_30[2], Turb_zo_30[2]) write.csv(Perc_zo_30, "Perc_zo_30.csv") write.csv(ENVI_zo_30, "ENVI_zo_30.csv") # winter Min_30 = cut(Chla_wi$Datetime, breaks = "30 min") chla_wi_30 = aggregate(Chl.a.ug.l. ~ Min_30, Chla_wi, mean) Temp_wi_30 = aggregate(Temp..deg.C. ~ Min_30, Chla_wi, mean) Turb_wi_30 = aggregate(Turb...M.FTU. ~ Min_30, Chla_wi, mean) Min_30 = cut(Flow_wi$Datetime, breaks = "30 min") Flow_wi_30 = aggregate(Speed ~ Min_30, Flow_wi, mean) Min_30 = cut(Perc_wi$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "30 min") Perc1_wi = aggregate(Perc.1 ~ Min_30, Perc_wi, mean) Perc6_wi = aggregate(Perc.6 ~ Min_30, Perc_wi, mean) Perc7_wi = aggregate(Perc.7 ~ Min_30, Perc_wi, mean) Perc8_wi = aggregate(Perc.8 ~ Min_30, Perc_wi, mean) Perc_wi_30 = cbind(Perc1_wi, Perc6_wi[2], Perc7_wi[2], Perc8_wi[2]) ENVI_wi_30 = cbind(chla_wi_30, Temp_wi_30[2], Turb_wi_30[2]) write.csv(Perc_wi_30, "Perc_wi_30.csv") write.csv(ENVI_wi_30, "ENVI_wi_30.csv") # 30-minute averaged: tide vs. flow (regression curves) plot(Flow_wi_30$Speed ~ Autumn_tide_30$Tide, ann = FALSE, las = 1) mtext(side = 3, "Autumn - Voordelta") mtext(side = 2, line = 3, expression("Flow velocity (m s"^{-1}*")")) mtext(side = 1, line = 3, "Tidal level (m)") model1 <-lm(Flow_wi_30$Speed ~ Autumn_tide_30$Tide+I(Autumn_tide_30$Tide^2)) plotdata <- cbind(Autumn_tide_30$Tide, predict(model1)) lines(plotdata[order(Autumn_tide_30$Tide),], col = "blue", lwd = 1.5) mtext(paste("R2 =",round(summary(model1)$r.squared,3)), side=3, line=-1.5, adj=0.05, cex = 1) mtext(paste("beta =",round(model1$coefficients[2],3)), side=3, line=-2.5, adj=0.05, cex = 1) mtext(paste("p =",round(summary(model1)$coef[2,4],5)), side=3, line=-3.5, adj=0.05, cex = 1) # Betonbak Min_30 = cut(Chla_BB$Datetime, breaks = "30 min") chla_BB_30 = aggregate(Chl.a.ug.l. ~ Min_30, Chla_BB, mean) Temp_BB_30 = aggregate(Temp..deg.C. ~ Min_30, Chla_BB, mean) Turb_BB_30 = aggregate(Turb...M.FTU. ~ Min_30, Chla_BB, mean) Min_30 = cut(Beton_tide$Datetime, breaks = "30 min") # The tide column is not a numeric # convert first with as.character(), then as.numeric() Beton_tide$Tide = as.numeric(as.character(Beton_tide$Tide)) Beton_tide_BB_30 = aggregate(Tide ~ Min_30, Beton_tide, mean) Beton_tide_BB_30$Tide = Beton_tide_BB_30$Tide/100 Min_30 = cut(Perc_BB$Datetime,format='%Y-%M-%d %H:%M:%S',tz="CET", breaks = "30 min") Perc1_BB = aggregate(Perc.1 ~ Min_30, Perc_BB, mean) Perc2_BB = aggregate(Perc.2 ~ Min_30, Perc_BB, mean) Perc3_BB = aggregate(Perc.3 ~ Min_30, Perc_BB, mean) Perc4_BB = aggregate(Perc.4 ~ Min_30, Perc_BB, mean) Perc5_BB = aggregate(Perc.5 ~ Min_30, Perc_BB, mean) Perc6_BB = aggregate(Perc.6 ~ Min_30, Perc_BB, mean) Perc7_BB = aggregate(Perc.7 ~ Min_30, Perc_BB, mean) Perc8_BB = aggregate(Perc.8 ~ Min_30, Perc_BB, mean) Perc_BB_30 = cbind(Perc1_BB, Perc2_BB[2], Perc3_BB[2], Perc4_BB[2], Perc5_BB[2], Perc6_BB[2], Perc7_BB[2], Perc8_BB[2]) ENVI_BB_30 = cbind(chla_BB_30, Temp_BB_30[2], Turb_BB_30[2]) # check if 30-minute rolling average vs. 30-minute averaged makes a difference (it doesn't really) ################################################################################# # This is the 30-minute MOVING average comparisons ################################################################################# par(mfrow = c(2,1), mar = c(2,6,2,6)) plot(ma(Perc_zo$Perc.5, order = 180, centre = TRUE), type = "l", las = 1, ylab = "", ylim = c(0, 1)) mtext(side = 2, line = 3, paste("Gape", i, sep=" ")) par(new=TRUE) plot(ma(Chla_zo$Chl.a.ug.l., order = 30, centre = TRUE), type = "l", col = "red", axes = F, xlab = NA, ylab = NA, lty = 3) mtext(side = 4, line = 3, expression( mu~"g L"^{-1})) axis(side = 4, at = seq(from = min(ENVI_zo_30$Chl.a.ug.l.), to = max(ENVI_zo_30$Chl.a.ug.l.), length.out = 5), labels = round(seq(from = min(ENVI_zo_30$Chl.a.ug.l.), to = max(ENVI_zo_30$Chl.a.ug.l.), length.out = 5), 0), las = 1) plot(Perc_zo_30$Perc.5, type = "l", las = 1, ylab = "", ylim = c(0, 1)) mtext(side = 2, line = 3, paste("Gape", i, sep=" ")) par(new=TRUE) plot(ENVI_zo_30$Chl.a.ug.l., type = "l", col = "red", axes = F, xlab = NA, ylab = NA, lty = 3) mtext(side = 4, line = 3, expression( mu~"g L"^{-1})) axis(side = 4, at = seq(from = min(ENVI_zo_30$Chl.a.ug.l.), to = max(ENVI_zo_30$Chl.a.ug.l.), length.out = 5), labels = round(seq(from = min(ENVI_zo_30$Chl.a.ug.l.), to = max(ENVI_zo_30$Chl.a.ug.l.), length.out = 5), 0), las = 1) ################################################################################# ################################################################################# # ------------------------------------ #### Betonbak vs. tidal # ------------------------------------ # all oysters were placed into betonbak 3 (1 closest to Delta Hall; 12 NIOZ meso) # The NAP height of the top of betonbak 3 is: 0.087 m # Raw data of gaping is too noisy; use 30-minute averaged or do a ma() Beton_tide$Tide = as.numeric(as.character(Beton_tide$Tide)) Beton_tide$Tide_m = Beton_tide$Tide / 100 par(mfrow = c(4,1), mar = c(2,6,2,6)) for(i in 2:9){ plot(Perc_BB_30[ ,i], type = "l", las = 1, ylab = "", ylim = c(-1, 1)) mtext(side = 2, line = 3, paste(names(Perc_BB_30)[i])) par(new=TRUE) plot(Beton_tide$Tide_m, ylim = c(-2, 4), axes = F, xlab = NA, ylab = NA, type = "l", col = "red", lty = 3) mtext(side = 4, line = 3, "m") axis(side = 4, at = seq(from = min(Beton_tide$Tide_m), to = max(Beton_tide$Tide_m), length.out = 5), labels = round(seq(from = min(Beton_tide$Tide_m), to = max(Beton_tide$Tide_m), length.out = 5), 0), las = 1) abline(h = 0.87, lty = 2, col = "blue", lwd = 1.5) } legend("topleft", c("Gape", "Tide"),bty="n",cex=1.25, lty = c(1, 3), col = c("black", "red"), lwd = c(2,2)) Position(function(x) x > 0.87, Beton_tide$Tide_m) highTide = Beton_tide$Tide_m > 0.87 #correct for the polygon (last point needs to match value of first point) polyTide = highTide * 100000 - 100000/2 polyTide = c(polyTide, polyTide[1]) tideSwitchPoints = c(1, cumsum(rle(highTide)$lengths) + 1) tideSwitchPoints = tideSwitchPoints[-length(tideSwitchPoints)] tideState = rle(highTide)$values switchTimes = Beton_tide$Datetime[tideSwitchPoints] findNearestTideState = function(x, switchTimes, tideState){ diff = x - switchTimes diff[diff > 0] = NA return(tideSwitchValues[which.min(abs(diff))]) } Perc_BB_30$Min_30 = as.POSIXct(Perc_BB_30$Min_30) Perc_BB_30$tideState = as.logical(sapply(Perc_BB_30$Min_30, FUN = findNearestTideState, switchTimes = switchTimes, tideState = tideState)) Perc_BB_30$tideState = !Perc_BB_30$tideState plot(Perc_BB_30$Perc.1[Perc_BB_30$tideState] ~ Perc_BB_30$Min_30[Perc_BB_30$tideState], pch = 22, bg = 4) points(Perc_BB_30$Perc.1[!Perc_BB_30$tideState] ~ Perc_BB_30$Min_30[!Perc_BB_30$tideState], pch = 21) Perc_BB_30$tideState[is.na(Perc_BB_30$tideState)] = T tideState = rle(Perc_BB_30$tideState)$value tideStarts = c(1, cumsum(rle(Perc_BB_30$tideState)$lengths) + 1) tideStarts = tideStarts[-length(tideStarts)] tideStarts = tideStarts[tideState == T] switchTimes = Perc_BB_30$Min_30[tideStarts] x = Perc_BB_30$Min_30[3] findTimeSinceTideOnset = function(x, switchTimes){ diff = as.numeric(difftime(x, switchTimes, units = 'hours')) diff = abs(diff[diff >= 0]) return(min(diff)) } Perc_BB_30$timeSinceTideOnset.hours = sapply(Perc_BB_30$Min_30, FUN = findTimeSinceTideOnset, switchTimes = switchTimes) Perc_BB_30 = cbind(Perc_BB_30, Beton_tide_BB_30$Tide) names(Perc_BB_30) = c("Min_30", "Oyster 1", "Oyster 2", "Oyster 3", "Oyster 4", "Oyster 5", "Oyster 6", "Oyster 7", "Oyster 8", "timeSinceTideOnset.hours", "Beton_tide_BB_30.Tide") # save the df to a file to load and make figures: write.csv(Perc_BB_30, "Perc_BB_30.csv")