data-raw/avasoilmoisture.R

In rutujact/avasoilmoisture: Compiles soil moisture data from Matteo's BCI Ava Tower TDR probes.

##-------------------
## Author: Rutuja Chitra-Tarak
## Date: 10/20/2018
## Title: Ava soil moisture provided by Matteo compiled
##-------------------
rm(list = ls())
if (!require("pacman")) install.packages("pacman"); library(pacman)
pacman::p_load(tidyverse, scales, cowplot, devtools, usethis)

# graphics info
tex <- element_text(size = 16, face = "plain") # , family = "gara"
my.theme <-  theme(axis.text = tex, axis.title = tex,
                   title = tex, legend.title = tex, legend.text = tex, strip.text.y = tex, strip.text.x = tex)
my.bg <- theme_bw()
my.adjust <- theme(axis.title.y = element_text(vjust = 1), axis.title.x = element_text(vjust = -0.6), title = element_text(vjust = 2))

path1 = "data-raw/AVA_SOIL"
all.files1 <- list.files(
  path =
    path1,
  pattern = ".csv",
  recursive = FALSE,
  full.names = TRUE
)

## files 2012-2015 only have data for the three vertical probe that covers depth 0-15 cm.
## While files after that have data for the above three as wells as three horizontal probes

total1 = length(all.files1);
# create progress bar
pb1 <- txtProgressBar(min = 0, max = total1, style = 3)
for (i in 1:length(all.files1)) {
  new.file <- read.csv(
    all.files1[i],
    skip = 4, ## remove details
    na.strings = c("NA", ""),
    header = F,
    row.names = NULL,
    check.names = F
  )
  if(i %in% 1:4){
    new.mod <- new.file %>% select(V1, V4:V6) %>%
      mutate_all(as.character) %>%
      mutate_at(vars(V4:V6), as.numeric)
  } else {
    new.mod <- new.file %>% select(V1, V4:V9) %>%
      mutate_all(as.character) %>%
      mutate_at(vars(V4:V9), as.numeric)
  }
  if (i == 1) {
    ava <- new.mod
  } else {
    ava <- bind_rows(ava, new.mod)
  }
  Sys.sleep(0.1)
  # update progress bar
  setTxtProgressBar(pb1, i)
  ## add headers
}
#
# colnames(ava) <- c("date.time", c("Depth_15.1", "Depth_15.2", "Depth_15.3",
#                                   "Depth_10", "Depth_40", "Depth_100"))

## Matteo's soil moisture data from Ava tower. This one is by the north side of BCI plot
## per Matteo's email dated Sep 17, 2018
## what you need is TDR_Avg(4) , TDR_Avg(5), TDR_Avg(6) which corresponds to 10 cm, 40 cm amd 100 cm.
## "If I remember well what you are interested in is the saturation index (0-100%)
# not the volumetric soil moisture (m3/m3), so the easy thing would [be]
# just to divide each series by its maximum value in the whole record
# (which should be representative of saturation conditions).


# tdr_model(x) = p1*x^3 + p2*x^2 + p3*x + p4
# Coefficients (with 95% confidence bounds):
# p1 =    0.001375  (0.0002419, 0.002508)
# p2 =     -0.1161  (-0.2152, -0.01708)
# p3 =       3.277  (0.4036, 6.151)
# p4 =      -30.57  (-58.23, -2.915)

head(ava); tail(ava)
sat.tdr <- ava
sat.tdr[ , -1] <- apply(ava[ ,-1], 2, function(x){x/max(x, na.rm = T)})
summary(ava)
summary(sat.tdr)

tdr <- gather(ava, key = depth, value = tdr, -V1) %>%
  mutate(date.time = lubridate::mdy_hm(V1)) %>%
  mutate(date = as.Date(date.time)) %>%
  select(-V1) %>%
  mutate(depth = plyr::mapvalues(depth, c(paste0("V", 4:9)),
                                 c("0-15_1", "0-15_2", "0-15_3", "10", "40", "100"))) %>%
  mutate(tdr = as.numeric(tdr),
         swc = 0.001375*tdr^3 -0.1161*tdr^2 + 3.277*tdr -30.57)

tdr$depth.date.time <- paste(tdr$depth, tdr$date.time, sep = ".")
tdr <- tdr[!duplicated(tdr$depth.date.time),]


swc.forwide <- tdr %>%select(date.time, depth, swc)
swc.wide <- spread(tdr.forwide, key = depth, value = swc)
sat.swc <- swc.wide
sat.swc[, -1] <- apply(swc.wide[ ,-1], 2, function(x){x/max(x, na.rm = T)})

sat.index <- gather(sat.swc, key = depth, value = sat.index, -date.time) %>%
  mutate(depth.date.time = paste(depth, date.time, sep = ".")) %>%
  select(-depth, -date.time)

tdr <- tdr %>% full_join(sat.index, by = "depth.date.time") %>%
  select(-depth.date.time) %>%
  select(date.time, date, depth, tdr, swc, sat.index)
head(tdr)

tdr.daily <- tdr %>% select(-date.time) %>%
  group_by(date, depth) %>%
  summarise_all(list(~mean(.)))

head(tdr)

horizontal <- tdr %>% subset(depth %in% c(10, 40, 100)) %>%
  subset(!is.na(swc)) %>%
  mutate(depth = as.numeric(depth))
vertical <- tdr %>% subset(!depth %in% c(10, 40, 100)) %>%
  subset(!is.na(swc))

write.csv(tdr,
          file.path("data-raw/tdr.csv"),
          row.names = FALSE)
write.csv(tdr.daily,
          file.path("data-raw/tdr.csv"),
          row.names = FALSE)
usethis::use_data(tdr, overwrite = TRUE)
usethis::use_data(tdr.daily, overwrite = TRUE)

usethis::use_data(horizontal, overwrite = TRUE)
usethis::use_data(vertical, overwrite = TRUE)

devtools::document()
# SWC
p1 <- ggplot(tdr, aes(x = date.time, y = swc, colour = depth)) +
  geom_line(aes(group = depth)) +
  ylab("Soil Moisture (m3/m3)") + xlab("Date") + #ylim(0.12, 0.8) +
  scale_x_datetime(date_breaks = "6 month", labels = function(x) format(x, "%d/%m/%y")) +
  my.theme + my.bg + my.adjust + theme(panel.grid.major.x = element_blank()) +
  theme(legend.position = c(0.95, 0.85), legend.background = element_rect(fill = "transparent")) +
  labs(colour = "Depth (cm)") +
  ggtitle("BCI: Derived Soil Water Content")
p1
ggsave(file.path("figures/Derived Soil Water Content all sensors.jpeg"), height = 7, width = 12, units='in')
p1 %+% subset(tdr, date > as.Date("2016-01-01"))
ggsave(file.path("figures/Derived Soil Water Content all sensors_common_period.jpeg"), height = 7, width = 12, units='in')

p1 %+% horizontal + aes(colour = as.factor(depth)) +
  theme(legend.position = c(0.9, 0.85)) + labs(colour = "Depth (cm)") +
  ggtitle("BCI: Derived SWC for depths 10, 40 and 100 cm")
ggsave(file.path("figures/Derived SWC for depths 10, 40 and 100 cm.jpeg"), height = 7, width = 12, units='in')

p1 %+% vertical + aes(colour = as.factor(depth)) +
  theme(legend.position = c(0.95, 0.9)) + labs(colour = "Depth (cm)") +
  ggtitle("BCI: Derived SWC for three vertical probes covering 0-15 cm")
ggsave(file.path("figures/Derived SWC for three vertical probes covering 0-15 cm.jpeg"), height = 7, width = 12, units='in')

## Saturation Index
p1 %+% horizontal + aes(x= date.time, y = sat.index, colour = as.factor(depth)) +
  ylab("Saturation Index (Unitless)") + xlab("Date") + #ylim(0.12, 0.8) +
  theme(legend.position = c(0.9, 0.85)) + labs(colour = "Depth (cm)") +
  ggtitle("BCI: Saturation Index for depths 10, 40 and 100 cm")
ggsave(file.path("figures/Saturation Index for depths 10, 40 and 100 cm.jpeg"), height = 7, width = 12, units='in')

## TDR
p1 %+% horizontal + aes(x= date.time, y = tdr, colour = as.factor(depth)) +
  ylab("TDR") + xlab("Date") + #ylim(0.12, 0.8) +
  theme(legend.position = c(0.9, 0.85)) + labs(colour = "Depth (cm)") +
  ggtitle("BCI: TDR data for depths 10, 40 and 100 cm")
ggsave(file.path("figures/TDR data for depths 10, 40 and 100 cm.jpeg"), height = 7, width = 12, units='in')