# Code by Benjamin Michael Marshall and Cameron Wesley Hodges 2021-08-01 for publication of: 
# title = {Malayan kraits (Bungarus candidus) show affinity to anthropogenic structures in a human dominated landscape}
## authors= {Cameron Wesley Hodges, Benjamin Michael Marshall, Jacques George Hill III, and Colin Thomas Strine}
## year = {2021}

activeData <- read.csv("Shelter_camera_trapping_activity_data.csv")

unique(activeData$Behavior)

library(ggplot2)
library(dplyr)
library(scales)
library(lubridate)
library(scico)

activeData$timestamp <- as.POSIXct(activeData$Scan.time,
                                   "%H:%M",  tz = "Asia/Bangkok")

library(hms)
activeData$time <- hms::hms(second(activeData$timestamp),
                            minute(activeData$timestamp),
                            hour(activeData$timestamp))
activeData$time <- as.POSIXct(activeData$time)


#     "university class period hours" bar with geom_segment (08:00-16:30)
#     "standard restaurant open hours" bar with geom_segment (08:00-21:00)
#     "bar open hours" bar with geom_segment (18:00-02:00)

annotationLayer <- 
  data.frame(
    label = c("University classes",
              "Restaurants open",
              "Bars open",
              "Bars open"
    ),
    start = as.POSIXct(c("08:00",
                         "08:00",
                         "18:00",
                         NA), format = "%H:%M",  tz = "Asia/Bangkok"),
    end = as.POSIXct(c("16:30",
                       "21:00",
                       NA,
                       "02:00"),
                     format = "%H:%M",  tz = "Asia/Bangkok"),
    ymin = c(-0.1, -0.2, -0.3, -0.30),
    ymax = c(-0.18, -0.28, -0.38, -0.38)
  )
annotationLayer$start <- hms::hms(second(annotationLayer$start),
                                  minute(annotationLayer$start),
                                  hour(annotationLayer$start))
annotationLayer$start <- as.POSIXct(annotationLayer$start)
annotationLayer$end <- hms::hms(second(annotationLayer$end),
                                minute(annotationLayer$end),
                                hour(annotationLayer$end))
annotationLayer$end <- as.POSIXct(annotationLayer$end)

annotationLayer[3,3] <- Inf
annotationLayer[4,2] <- -Inf


# 05:40 to 06:40
# 18:45 to 17:50

minDate1 <- as.POSIXct(c("1970-01-01 05:40"), format = "%Y-%m-%d %H:%M")
maxDate1 <- as.POSIXct(c("1970-01-01 06:40"), format = "%Y-%m-%d %H:%M")
dates1 <- seq(minDate1, maxDate1, by = "min")

minDate2 <- as.POSIXct(c("1970-01-01 17:50"), format = "%Y-%m-%d %H:%M")
maxDate2 <- as.POSIXct(c("1970-01-01 18:45"), format = "%Y-%m-%d %H:%M")
dates2 <- seq(minDate2, maxDate2, by = "min")

lefts <- seq(minDate1, maxDate1, 0.5)
rights <- c(lefts[-1], maxDate1 + 0.5)
maxs <- rep(Inf, length(lefts))
mins <- rep(-Inf, length(lefts))
fill <- seq(0, 1, length.out = length(lefts))
shadeDf <- data.frame(lefts, rights, maxs, mins, fill)

lefts <- seq(minDate2, maxDate2, 0.5)
rights <- c(lefts[-1], maxDate2 + 0.5)
maxs <- rep(Inf, length(lefts))
mins <- rep(-Inf, length(lefts))
fill <- seq(0, 1, length.out = length(lefts))
shadeDf2 <- data.frame(lefts, rights, maxs, mins, fill = sort(fill, decreasing = TRUE))

shadeAll <- rbind(shadeDf, shadeDf2)

# fill in middle
lefts <- seq(maxDate1, minDate2, by = "min")
rights <- c(lefts[-1], minDate2 + 0.5)
maxs <- rep(Inf, length(lefts))
mins <- rep(-Inf, length(lefts))
fill <- rep(1, length(lefts))
shadeDfMiddle <- data.frame(lefts, rights, maxs, mins, fill)

shadeAll <- rbind(shadeAll, shadeDfMiddle)

bucaActColour <- scico(1, palette = "roma", direction = -1)

activeData %>% 
  filter(!Behavior %in% c("Sheltering", "Unknown", "")) %>% 
  ggplot() + 
  geom_rect(data = shadeAll,
            aes(xmin = lefts, xmax = rights, ymin = mins, ymax = maxs,
                fill = fill, alpha = fill)) +
  geom_density(aes(x = time, y = ..scaled..),
               fill = bucaActColour,
               colour = bucaActColour,
               alpha = 0.75) +
  geom_point(aes(x = time, y = -0.05), pch = "|", size = 4, alpha = 0.5,
             colour = bucaActColour) +
  # geom_segment(data = annotationLayer,
  #              aes(x = start, xend = end, y = y, yend = y)) +
  geom_rect(data = annotationLayer,
            aes(xmin = start, xmax = end, ymin = ymin, ymax = ymax)) +
  geom_text(data = annotationLayer,
            aes(x = start, y = (ymin+ymax)/2, label = label),
            colour = "white", hjust = -0.02, vjust = 0.5, fontface = 2,
            size = 3) +
  annotate("text", x = as.POSIXct("1970-01-01 12:15:00"), y = 0.95,
           label = "DAYLIGHT", size = 6, fontface = 2,
           colour = "black", hjust = 0.5, alpha = 0.65) +
  annotate("text", x = as.POSIXct("1970-01-01 00:20:00"), y = 0.03,
           label = "BUCA activity", size = 4, fontface = 2,
           colour = bucaActColour, hjust = 0, vjust = 0, alpha = 0.75) +
  labs(x = "Time (24 hour)") +
  scale_fill_scico(palette = "nuuk", begin = 0.6) +
  scale_alpha_continuous(range = c(0.1, 0.5)) +
  scale_x_datetime(breaks = date_breaks("1 hours"),
                   labels = date_format("%H"),
                   limits = as.POSIXct(c("1970-01-01 00:00:01",
                                         "1970-01-01 23:59:59"),
                                       format = "%Y-%m-%d %H:%M:%S"),
                   expand = c(0,0)) +
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(),
        axis.title.x = element_text(face = 2),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank())

ggsave("./Figures/Activity Plot.png", units = "mm", dpi = 300,
       width = 160, height = 90)
ggsave("./Figures/Activity Plot.pdf", units = "mm",
       width = 160, height = 90)
