---
title: "Consolidate Data"
author: "Sydney Collins"
date: "2023-06-10"
output: html_document
editor_options: 
  chunk_output_type: console
---

#Load Packages
```{r Libraries}
rm(list=ls())
options(scipen = 999)

library(ncdf4)
library(tidyverse)
library(lubridate)
library(data.table)
library(chillR)
library(moonlit)
library(suncalc)
```

#Cloud Cover Dataset

Data from https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-cloud-properties?tab=form

```{r Cloud Cover Dataset}
#Formatting the cloud cover data
file_thesis <- list.files("data/CloudCover", pattern = "*.nc", all.files = TRUE, full.names = TRUE)
cloud <- data.frame()

for(i in 1:length(file_thesis)){
  cdf <- nc_open(file_thesis[i])
  a <- as.data.frame(ncvar_get(cdf, varid="cfc"))
  x <- a[29,25]
  y <- file_thesis[i]
  z <- cbind(x, y)
  cloud <- rbind(cloud, z)
}

cloud$y <- gsub("000040019AVPOS01GL_area_subset.nc", "", cloud$y)
cloud$y <- gsub("data/CloudCover/5be64465-46fb-41dd-9a7d-4dc86e088f3c-CFCdm", "", cloud$y)
cloud$y <- gsub("000000219AVPOSE1GL_area_subset.nc", "", cloud$y)
cloud$y <- gsub("data/CloudCover/2929ea4c-02bc-4441-9ce6-dc3805119109-CFCdm", "", cloud$y)
cloud$y <- gsub("data/CloudCover/CFCdm", "", cloud$y)
cloud$y <- gsub("000040019AVPOS01GL.nc", "", cloud$y)

names(cloud) <- c("CloudCover", "Date")
cloud$Date <- gsub("^(.{4})(.*)$", "\\1-\\2", cloud$Date)
cloud$Date <- gsub("^(.{7})(.*)$", "\\1-\\2", cloud$Date)
cloud$Date <- as.Date(cloud$Date)
cloud$CloudCover <- as.integer(cloud$CloudCover)

save(cloud, file = "data/cloud.Rda")
```

#Moon Dataset
```{r Moon Phase}
#Average moon phase for each night throughout the fledging season
start_date17 <- as.POSIXct("2017-09-13 01:00:00")
start_date18 <- as.POSIXct("2018-09-13 01:00:00")
start_date21 <- as.POSIXct("2021-09-13 01:00:00")
start_date22 <- as.POSIXct("2022-09-13 01:00:00")

dates17 <- as.data.frame(seq(start_date17, by = "day", length.out = 62)) %>%
  rename(DateTime = `seq(start_date17, by = "day", length.out = 62)`)
dates17 <- calculateMoonlightStatistics(lat = 47.26265, lon = -52.77187, date = dates17$DateTime, e = 0.26, t = "1 hour", timezone = "America/St_Johns")

dates18 <- as.data.frame(seq(start_date18, by = "day", length.out = 62)) %>%
  rename(DateTime = `seq(start_date18, by = "day", length.out = 62)`)
dates18 <- calculateMoonlightStatistics(lat = 47.26265, lon = -52.77187, date = dates18$DateTime, e = 0.26, t = "1 hour", timezone = "America/St_Johns")

dates21 <- as.data.frame(seq(start_date21, by = "day", length.out = 62)) %>%
  rename(DateTime = `seq(start_date21, by = "day", length.out = 62)`)
dates21 <- calculateMoonlightStatistics(lat = 47.26265, lon = -52.77187, date = dates21$DateTime, e = 0.26, t = "1 hour", timezone = "America/St_Johns")

dates22 <- as.data.frame(seq(start_date22, by = "day", length.out = 62)) %>%
  rename(DateTime = `seq(start_date22, by = "day", length.out = 62)`)
dates22 <- calculateMoonlightStatistics(lat = 47.26265, lon = -52.77187, date = dates22$DateTime, e = 0.26, t = "1 hour", timezone = "America/St_Johns")

phasedates <- rbind(dates17, dates18, dates21, dates22) %>%
  mutate(Phase = as.factor(ifelse(meanMoonPhase < 0.25, "first quarter", ifelse(meanMoonPhase < 0.5, "second quarter", ifelse(meanMoonPhase < 0.75, "third quarter", "fourth quarter"))))) %>%
  mutate(phasecat = sjmisc::group_var(meanMoonPhase, size = 0.05, as.num = FALSE)) %>%
  mutate(phasecat = as.numeric(phasecat)*5) %>%
  mutate(Date = as.Date(date))

save(phasedates, file = "data/phasedates.Rda")
```

```{r Moon Illumination}
#Incident moon illumination for each night throughout the fledging season
start_datetime17 <- as.POSIXct("2017-09-13 00:00")
start_datetime18 <- as.POSIXct("2018-09-13 00:00")
start_datetime21 <- as.POSIXct("2021-09-13 00:00")
start_datetime22 <- as.POSIXct("2022-09-13 00:00")
`%ni%` <- Negate(`%in%`)

datetime17 <- as.data.frame(seq(start_datetime17, by = "min", length.out = 89340)) %>%
  rename(DateTime = `seq(start_datetime17, by = "min", length.out = 89340)`) %>%
  mutate(hour = hour(DateTime))
datetime17 <- calculateMoonlightIntensity(lat = 47.26265, lon = -52.77187, date = datetime17$DateTime, e = 0.26) %>%
  filter(night == "TRUE")

datetime18 <- as.data.frame(seq(start_datetime18, by = "min", length.out = 89340)) %>%
  rename(DateTime = `seq(start_datetime18, by = "min", length.out = 89340)`) %>%
  mutate(hour = hour(DateTime))
datetime18 <- calculateMoonlightIntensity(lat = 47.26265, lon = -52.77187, date = datetime18$DateTime, e = 0.26) %>%
  filter(night == "TRUE")

datetime21 <- as.data.frame(seq(start_datetime21, by = "min", length.out = 89340)) %>%
  rename(DateTime = `seq(start_datetime21, by = "min", length.out = 89340)`) %>%
  mutate(hour = hour(DateTime))
datetime21 <- calculateMoonlightIntensity(lat = 47.26265, lon = -52.77187, date = datetime21$DateTime, e = 0.26) %>%
  filter(night == TRUE)

datetime22 <- as.data.frame(seq(start_datetime22, by = "min", length.out = 89340)) %>%
  rename(DateTime = `seq(start_datetime22, by = "min", length.out = 89340)`) %>%
  mutate(hour = hour(DateTime))
datetime22 <- calculateMoonlightIntensity(lat = 47.26265, lon = -52.77187, date = datetime22$DateTime, e = 0.26) %>%
  filter(night == "TRUE")

moondates <- rbind(datetime17, datetime18, datetime21, datetime22) %>%
  filter(night == "TRUE") %>%
  mutate(MoonUp = as.factor(ifelse(moonAltDegrees > 0, "UP", "DOWN"))) %>%
  mutate(Phase = as.factor(ifelse(moonPhase < 0.25, "first quarter", ifelse(moonPhase < 0.5, "second quarter", ifelse(moonPhase < 0.75, "third quarter", "fourth quarter"))))) %>%
  mutate(phasecat = sjmisc::group_var(moonPhase, size = 0.05, as.num = FALSE)) %>%
  mutate(phasecat = as.numeric(phasecat)*5) %>%
  mutate(illumcat = sjmisc::group_var(moonlightModel, size = 0.05, as.num = FALSE)) %>%
  mutate(illumcat = ifelse(moonlightModel == 0, "0", illumcat)) %>%
  mutate(illumcat = as.numeric(illumcat)*5)

save(moondates, file = "data/moondates.Rda")
```

#Chick Fledging Dataset
```{r Chick Fledge Age}
#Load data
fledgeDate <- read.csv("data/all_chicks.csv") %>%
  arrange(BandNo, DateTime) %>%
  group_by(BandNo) %>%
  filter(n() > 1) %>%
  slice_tail() %>%
  ungroup() %>%
  select(year, Plot, BurrowID, BandNo, DateTime) %>%
  rename(Year = year)

#Estimate age at last read using unpublished formula derived by R.A. Mauck
wl <- read.csv("data/chick_banding.csv") %>%
  filter(Tag.deployed != "") %>%
  select(Year, Date, Plot, BurrowID, BandNo, BirdWt, Wing) %>%
  left_join(fledgeDate, by = c("Year", "BandNo")) %>%
  mutate(Date = as.Date(Date), Fledge = as.Date(DateTime)) %>%
  mutate(DaysSinceBanding = as.numeric(Fledge - Date)) %>%
  mutate(AgeatFledge = 0.297*Wing + 13.521 + DaysSinceBanding) %>%
  mutate(FledgeWLmin = Wing + (DaysSinceBanding*(2.24-0.08)), FledgeWLmax = Wing +(DaysSinceBanding*(2.24+0.08))) %>%
  select(-Plot.y) %>%
  select(-BurrowID.y) %>%
  rename(Plot = Plot.x, BurrowID = BurrowID.x)

boxplot(wl$AgeatFledge)

#Eliminate birds that are not old enough to fledge (min 56 days)
nofledge <- wl %>% filter(AgeatFledge < 56)
nofledge <- as.character(nofledge$BandNo)
save(nofledge, file = "data/nofledge.Rda")
```

```{r Chick Fledging Dates}
##Format the data##
chick <- fledgeDate %>%
  mutate(DateTime = as.POSIXct(DateTime, format = "%Y-%m-%d %H:%M", tz = "America/St_Johns")) %>%
  mutate(yDate = yday(as.POSIXct(DateTime))) %>%
  mutate(Time = as.ITime(DateTime)) %>%
  mutate(ITime = as.numeric(Time)) %>%
  mutate(Date = as.Date(DateTime)) %>%
  filter(!(BandNo %in% nofledge ))

#scale the time to 6pm
chick_2 <- chick %>% filter(ITime < 60000)
chick <- chick %>% filter(ITime > 60000)
chick_2$ITime <- chick_2$ITime + 86400
chick <- rbind(chick, chick_2)
chick$ITime <- (chick$ITime - 64800)/60
chick$Year <- as.factor(chick$Year)
chick$Plot <- as.factor(chick$Plot)
```

#Combine Datasets
```{r Moon Fledge}
#Add the cloud cover data
chick <- left_join(chick, cloud, by = "Date") %>%
  arrange(DateTime)

#Add the sunset data
sunset1 <- as.data.frame(getSunlightTimes(date = chick$Date, lat = 47.26265,lon = -52.77187, tz = "America/St_Johns")) %>%
  mutate(Sunset = as.numeric(as.ITime(sunset))/(60*60))
chick$Sunset <- as.integer(sunset1$Sunset*60 - (18*60))
chick$TimeAfterSunset <- (chick$ITime - chick$Sunset)/60

#Moon data from moonlit
chick_moon <- calculateMoonlightIntensity(lat = 47.26265, lon = -52.77187, date = chick$DateTime, e = 0.26) %>%
  arrange(date)

chick_phase <- calculateMoonlightStatistics(lat = 47.26265, lon = -52.77187, date = chick$DateTime, e = 0.26, t = "1 hour", timezone = "America/St_Johns") %>%
  arrange(date)

chick <- cbind(chick, chick_moon[,4:10], chick_phase$meanMoonPhase) %>%
  rename(meanMoonPhase = 'chick_phase$meanMoonPhase') %>%
  mutate(MoonUp = as.factor(ifelse(moonAltDegrees > 0, "UP", "DOWN"))) %>%
  mutate(phasecat = sjmisc::group_var(meanMoonPhase, size = 0.05, as.num = FALSE)) %>%
  mutate(phasecat = as.numeric(phasecat)*5) %>%
  mutate(phasecat = ifelse(phasecat > 25, phasecat + 5, phasecat)) %>%
  mutate(illumcat = sjmisc::group_var(moonlightModel, size = 0.05, as.num = FALSE)) %>%
  mutate(illumcat = ifelse(moonlightModel < 0.00001, "0", illumcat)) %>%
  mutate(illumcat = as.numeric(illumcat)*5) %>%
  mutate(illumcat = ifelse(illumcat > 20, illumcat + 5, illumcat)) %>%
  mutate(Phase = as.factor(ifelse(meanMoonPhase < 0.25, "first quarter", ifelse(meanMoonPhase < 0.5, "second quarter", ifelse(meanMoonPhase < 0.75, "third quarter", "fourth quarter"))))) %>%
  left_join(wl %>% select(BandNo, AgeatFledge), by = "BandNo")

#Save the dataframes
save(chick, file = "data/chick.Rda")
```

#Strandings Data

Data from Burt, T. V. 2022, August 1. Influence of anthropogenic light on the attraction and mortality of Leach’s Storm-Petrels: Where, when, why, and which birds are stranding. Undergraduate dissertation, Memorial University of Newfoundland and Labrador, St. John’s.

```{r Strandings Data}
#Create and clean up the strandings dataset
#Want to get rid of days that had the lights on because this inflates the number of stranded birds
strand_BdV <- read.csv("data/Fishplant_2022_Feb_27_2023_SG.csv") %>%
  mutate(Date = as.POSIXlt(Date, format = "%m-%d-%Y")) %>%
  mutate(Year = Date$year + 1900) %>%
  filter(Year > 2018) %>%
  filter(Database %like% "Sherry's Data") %>%
  mutate(Stranded = as.numeric(Total_LESP)) %>%
  mutate(Lunar = as.integer(round(Lunar))) %>%
  dplyr::select(-(Plant_Lights_on.off:Number.of.Boxes)) %>%
  dplyr::select(-(Total_LESP:Database)) %>%
  drop_na() %>%
  filter(Stranded < 100)

strand_rep <- strand_BdV[rep(row.names(strand_BdV), strand_BdV$Stranded), 1:5]

strand <- calculateMoonlightStatistics(lat = 48.09780, lon = -52.89831, date = as.POSIXct(strand_rep$Date), e = 0.26, t = "1 hour", timezone = "America/St_Johns") %>%
  mutate(phasecat = sjmisc::group_var(meanMoonPhase, size = 0.05, as.num = FALSE)) %>%
  mutate(phasecat = as.numeric(phasecat)*5)

save(strand, file = "data/strand.Rda")
```