Argos Telemetry Data

BWTE 2013-2017

Contact: jmh09r@my.fsu.edu

date() 

## [1] "Mon Mar 19 19:19:39 2018"

Options

library(knitr)
library(rgl)
opts_knit$set(verbose = TRUE)
knit_hooks$set(webgl = hook_webgl)

Set Directory

setwd("D:/Patuxent/BWTE_ARamey/BWTE_ARamey_2013-2017_Argos_data-20180307T004644Z-001/BWTE_ARamey_2013-2017_Argos_data/Argos_data")

Load Packages:

suppressMessages(library(dplyr))
suppressMessages(library(data.table))
suppressMessages(library(readr))
suppressMessages(library(stringr))
suppressMessages(library(lubridate))
suppressMessages(library(RcppRoll))

## Warning: package 'RcppRoll' was built under R version 3.4.4

suppressMessages(library(changepoint))

## Warning: package 'changepoint' was built under R version 3.4.4

suppressMessages(library(ecp))

## Warning: package 'ecp' was built under R version 3.4.4

suppressMessages(library(trend))

## Warning: package 'trend' was built under R version 3.4.4

suppressMessages(library(ggplot2))
suppressMessages(library(gridExtra))

Load Data

x = "BWTeal_2013-2017_Ramey_argos_tabular_calibrated_sensor.csv"

df1 = as.data.frame(
          read_delim(x, 
                     col_names = TRUE, 
                     delim = ",", 
                     quoted_na = FALSE, 
                     trim_ws = TRUE))

## Warning: Missing column names filled in: 'X12' [12]

## Parsed with column specification:
## cols(
##   Data_DOI = col_character(),
##   Species = col_character(),
##   Animal_ID = col_character(),
##   Alive = col_character(),
##   Temperature_C = col_double(),
##   Voltage = col_double(),
##   Activity = col_integer(),
##   Program = col_integer(),
##   PTT = col_integer(),
##   Satellite = col_character(),
##   Message_Datetime_GMT = col_datetime(format = ""),
##   X12 = col_character()
## )

PointChange Analysis

Perm = 500 #permutations

Animal.IDs = levels(as.factor(df1$Animal_ID))

DeadShed.df = data.frame(matrix(nrow = length(Animal.IDs), ncol=16))

names(DeadShed.df) = c("Animal", "Act_c", "Act_d", "Temp_c","Temp_d", "mTmp_c", "mTmp_d", 
                       "UnqAct_c", "UnqAct_d", "UnqActR_c", "UnqActR_d", "Sat_c","Sat_d", 
                       "SatR_c","SatR_d","V_ln")

DeadShed.df$Animal = Animal.IDs

Loop

for(i in 1:length(Animal.IDs)){
  
          df1a.tmp = df1 %>%
                      filter(Animal_ID == Animal.IDs[i]) %>%
                      mutate(ObsDate = ymd_hms(Message_Datetime_GMT),
                             ObsDateD = as.Date(Message_Datetime_GMT))
          
          df1a.tmp = arrange(df1a.tmp, ObsDate) %>%
                      mutate(IDx = 1:nrow(df1a.tmp),
                             Activity = as.numeric(Activity),
                             A.dif = Activity - lag(Activity),
                             Year = year(ObsDate),
                             Month = month(ObsDate),
                             Week = week(ObsDate),
                             Day = day(ObsDate),
                             Hour = hour(ObsDate),
                             Seconds = second(ObsDate),
                             Sec.d = as.numeric(ObsDate, units = "seconds"),
                             Dur = Sec.d - lag(Sec.d),
                             Min.d = Dur/60,
                             Hr.d = Min.d/60,
                             Day.d = Hr.d/24)
          
           Cutp = br.test(df1a.tmp$Activity, m = Perm)
           
           Cutp2 = bu.test(df1a.tmp$Activity, m = Perm)
           
           Cutp3 = lanzante.test(df1a.tmp$Activity, method = "wilcox.test")
           
           Cutp4 = pettitt.test(df1a.tmp$Activity)
           
           Cutp5 = snh.test(df1a.tmp$Activity, m = Perm)
           
           SCut = max(c(Cutp$estimate, Cutp2$estimate, Cutp3$estimate, 
                        Cutp4$estimate, Cutp5$estimate), na.rm=T)
           
           
           
           CutpB = br.test(df1a.tmp$Temperature_C, m = Perm)
           
           CutpB1 = br.test(df1a.tmp$Temperature_C, m = Perm)
           
           CutpB2 = bu.test(df1a.tmp$Temperature_C, m = Perm)
           
           CutpB3 = lanzante.test(df1a.tmp$Temperature_C, method = "wilcox.test")
           
           CutpB4 = pettitt.test(df1a.tmp$Temperature_C)
           
           CutpB5 = snh.test(df1a.tmp$Temperature_C, m = Perm)
           
           SCutB = max(c(CutpB$estimate, CutpB2$estimate, CutpB3$estimate, 
                         CutpB4$estimate, CutpB5$estimate), na.rm=T)
           
           
           df1a.tmp$rmT= rollapply(df1a.tmp$Temperature_C, 2, mean, partial=T, align = "right")
            
           CutpBm = br.test(df1a.tmp$rmT, m = Perm)
           
           CutpBm1 = br.test(df1a.tmp$Activity, m = Perm)
           
           CutpBm2 = bu.test(df1a.tmp$rmT, m = Perm)
           
           CutpBm3 = lanzante.test(df1a.tmp$rmT, method = "wilcox.test")
           
           CutpBm4 = pettitt.test(df1a.tmp$rmT)
           
           CutpBm5 = snh.test(df1a.tmp$rmT, m = Perm)
           
           SCutBm = max(c(CutpBm$estimate, CutpBm2$estimate, CutpBm3$estimate, 
                          CutpBm4$estimate, CutpBm5$estimate), na.rm=T)
           
           
           
           Act.df = as.data.frame(df1a.tmp %>%
                     group_by(ObsDateD) %>%
                     summarise(Act = length(unique(Activity)),
                               Idx = max(IDx)))
           
           CutpC = br.test(Act.df$Act, m = Perm)
           
           CutpC2 = bu.test(Act.df$Act, m = Perm)
           
           CutpC3 = lanzante.test(Act.df$Act, method = "wilcox.test")
           
           CutpC4 = pettitt.test(Act.df$Act)
           
           CutpC5 = snh.test(Act.df$Act, m = Perm)
           
           SCutC = max(c(CutpC$estimate, CutpC2$estimate, CutpC3$estimate, 
                         CutpC4$estimate, CutpC5$estimate))
           
           GetCutUnqA = Act.df$Idx[SCutC]
           
           
           rl = rle(Act.df$Act)
           Act.df$repx = rep( rl$lengths != 1 , times = rl$lengths )
           Act.df$repx2 = ifelse(Act.df$repx == TRUE, 1, 0)
           
           if(sum(Act.df$repx2) == 0){
             
             DeadShed.df[i, 10] = 0
             DeadShed.df[i, 11] = 0 }

              else{

                   CutpCx = br.test(Act.df$repx2, m = Perm)
                   
                   CutpCx2 = bu.test(Act.df$repx2, m = Perm)
                   
                   CutpCx3 = lanzante.test(Act.df$repx2, method = "wilcox.test")
                   
                   CutpCx4 = pettitt.test(Act.df$repx2)
                   
                   CutpCx5 = snh.test(Act.df$repx2, m = Perm)
                   
                   SCutCx = max(c(CutpCx$estimate, CutpCx2$estimate, CutpCx3$estimate, 
                                  CutpCx4$estimate, CutpCx5$estimate))
                   
                   GetCutUnqR = Act.df$Idx[SCutCx]
                   
                   uActR.D = Act.df$ObsDateD[SCutCx]
                   
                   DeadShed.df[i, 10] = GetCutUnqR
                   DeadShed.df[i, 11] = strftime(uActR.D,  format = "%Y-%m-%d")} 
           
           
           
           Sat.df = as.data.frame(df1a.tmp %>%
                     group_by(ObsDateD) %>%
                     summarise(Cnt = length(ObsDate),
                               Idx = max(IDx)))
           
           CutpS = br.test(Sat.df$Cnt, m = Perm)
           
           CutpS2 = bu.test(Sat.df$Cnt, m = Perm)
           
           CutpS3 = lanzante.test(Sat.df$Cnt, method = "wilcox.test")
           
           CutpS4 = pettitt.test(Sat.df$Cnt)
           
           CutpS5 = snh.test(Sat.df$Cnt, m = Perm)
           
           SCutS = max(c(CutpS$estimate, CutpS2$estimate, CutpS3$estimate, 
                         CutpS4$estimate, CutpS5$estimate))
           
           GetCutSat = Sat.df$Idx[SCutS]
           
           
           
           rl = rle(Sat.df$Cnt)
           Sat.df$repx = rep( rl$lengths != 1 , times = rl$lengths )
           Sat.df$repx2 = ifelse(Sat.df$repx == TRUE, 1, 0)
           
           if(sum(Sat.df$repx2) == 0){
             
                DeadShed.df[i, 14] = 0 
                DeadShed.df[i, 15] = 0}

              else{

                   CutpCSx = br.test(Sat.df$repx2, m = Perm)
                   
                   CutpCSx2 = bu.test(Sat.df$repx2, m = Perm)
                   
                   CutpCSx3 = lanzante.test(Sat.df$repx2, method = "wilcox.test")
                   
                   CutpCSx4 = pettitt.test(Sat.df$repx2)
                   
                   CutpCSx5 = snh.test(Sat.df$repx2, m = Perm)
                   
                   SCutCSx = max(c(CutpCSx$estimate, CutpCSx2$estimate, CutpCSx3$estimate, CutpCSx4$estimate, CutpCSx5$estimate))
                   
                   GetCutSatR = Sat.df$Idx[SCutCSx] 
                   
                   SatR.D = Sat.df$ObsDateD[SCutCSx]
                   
                   DeadShed.df[i, 14] = GetCutSatR 
                   DeadShed.df[i, 15] = strftime(SatR.D,  format = "%Y-%m-%d")}
           

           Act.D = df1a.tmp$ObsDate[SCut]
           Temp.D = df1a.tmp$ObsDate[SCutB]
           mTemp.D = df1a.tmp$ObsDate[SCutBm]
           uAct.D = Act.df$ObsDateD[SCutC]
           Sat.D = Sat.df$ObsDateD[SCutS]
        
           DeadShed.df[i, 2] = SCut
           DeadShed.df[i, 3] = strftime(Act.D,  format = "%Y-%m-%d %H:%M:%S")
           DeadShed.df[i, 4] = SCutB
           DeadShed.df[i, 5] = strftime(Temp.D,  format = "%Y-%m-%d %H:%M:%S")
           DeadShed.df[i, 6] = SCutBm
           DeadShed.df[i, 7] = strftime(mTemp.D,  format = "%Y-%m-%d %H:%M:%S")
           DeadShed.df[i, 8] = GetCutUnqA
           DeadShed.df[i, 9] = strftime(uAct.D,  format = "%Y-%m-%d")
           DeadShed.df[i, 12] = GetCutSat 
           DeadShed.df[i, 13] = strftime(Sat.D,  format = "%Y-%m-%d")
           DeadShed.df[i, 16] = dim(df1a.tmp)[1]
           
}

write.csv(DeadShed.df, "./DeadShed031918.csv")

Add Column with Dead/Live Designation

for(i in 1:length(Animal.IDs)){
  
          DeathPoint = (DeadShed.df %>%
                       filter(Animal == Animal.IDs[i]))[,18]
          
          Tmp.mrk = df1 %>%
                      filter(Animal_ID == Animal.IDs[i]) %>%
                      mutate(ObsDate = ymd_hms(Message_Datetime_GMT),
                             ObsDateD = as.Date(Message_Datetime_GMT),
                             Live2 = NA)
          
          Tmp.mrk = arrange(Tmp.mrk, ObsDate) %>%
                    mutate(Idx = 1:nrow(Tmp.mrk))
          
          Tmp.mrk$Alive2 = ifelse(Tmp.mrk$Idx <= DeathPoint, "Alive", "Dead")
          
              if(i == 1){mrkOut.df = Tmp.mrk}
                  else{mrkOut.df = rbind(mrkOut.df, Tmp.mrk)}
          
}

write.csv(mrkOut.df, "./mrkd_BWTeal_2013-2017_argos_031918.csv") #w/ markup

mrkOut2.df = mrkOut.df %>%
              filter(Alive2 == "Alive")

write.csv(mrkOut2.df, "./clean_BWTeal_2013-2017_argos_031918.csv") #Copy with dead removed
    
   
#COunt of records removed
dim(mrkOut.df)[1]-dim(mrkOut2.df)[1]

## [1] 3859

Diagnostic Plots

Showing “DeadPoint” on orginal dataset.

for(Setx in 1:length(Animal.IDs)){

Plt.set = df1 %>%
          filter(Animal_ID == Animal.IDs[Setx]) %>%
          mutate(ObsDate = ymd_hms(Message_Datetime_GMT),
                 ObsDateD = as.Date(Message_Datetime_GMT))

Plt.set = arrange(Plt.set, ObsDate) %>%
          mutate(IDx = 1:nrow(Plt.set),
                 Activity = as.numeric(Activity))

DeathPoint = (DeadShed.df %>%
              filter(Animal == Animal.IDs[Setx]))[,18]


P1 = ggplot(Plt.set, aes(IDx,  Temperature_C)) +  
       geom_line() +
       geom_vline(xintercept=DeathPoint, col = "darkred", lwd = 1.25, linetype = "dotted") +
       ylab("Temperature") + 
       xlab("Time Step") + 
       ylim(0, 80) +
       ggtitle(paste(Animal.IDs[Setx])) +
       theme_classic() +
       guides(col = guide_legend(ncol = 16)) +
       theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
              strip.text.x = element_text(size=12, face="bold"),
              axis.title.y = element_text(size=16, face="bold"),
              plot.title = element_text(hjust = 0.5),
              legend.position="bottom")
      


P2 = ggplot(Plt.set, aes(IDx,  Activity)) +  
       geom_line() +
       geom_vline(xintercept=DeathPoint, col = "darkred", lwd = 1.25, linetype = "dotted") +
       ylab("Activity") + 
       xlab("Time Step") + 
       ylim(0, 100) +
       theme_classic() +
       ggtitle(paste(Animal.IDs[Setx])) +
       guides(col = guide_legend(ncol = 16)) +
       theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
              strip.text.x = element_text(size=12, face="bold"),
              axis.title.y = element_text(size=16, face="bold"),
              plot.title = element_text(hjust = 0.5),
              legend.position="bottom")



Act.plt = as.data.frame(Plt.set %>%
                     group_by(ObsDateD) %>%
                     summarise(Act = length(unique(Activity)),
                               Idx = max(IDx)))

DeathDate = Plt.set$ObsDateD[DeathPoint]

P3 = ggplot(Act.plt, aes(ObsDateD,  Act)) +  
       geom_bar(colour="black", stat="identity") +
       geom_vline(xintercept=DeathDate, col = "darkred", lwd = 1.25, linetype = "dotted") +
       ylab("Activity/Day") + 
       xlab("Day") + 
       #ylim(0, 100) +
       theme_classic() +
       ggtitle(paste(Animal.IDs[Setx])) +
       guides(col = guide_legend(ncol = 16)) +
       theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
              strip.text.x = element_text(size=12, face="bold"),
              axis.title.y = element_text(size=16, face="bold"),
              plot.title = element_text(hjust = 0.5),
              legend.position="bottom")


 Sat.plt = as.data.frame(Plt.set %>%
                     group_by(ObsDateD) %>%
                     summarise(Cnt = length(ObsDate),
                               Idx = max(IDx)))
 
 P4 = ggplot(Sat.plt, aes(ObsDateD,  Cnt)) +  
       geom_bar(colour="black", stat="identity") +
       geom_vline(xintercept=DeathDate, col = "darkred", lwd = 1.25, linetype = "dotted") +
       ylab("Transmit/Day") + 
       xlab("Day") + 
       ylim(0, 150) +
       theme_classic() +
       ggtitle(paste(Animal.IDs[Setx])) +
       guides(col = guide_legend(ncol = 16)) +
       theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
              strip.text.x = element_text(size=12, face="bold"),
              axis.title.y = element_text(size=16, face="bold"),
              plot.title = element_text(hjust = 0.5),
              legend.position="bottom")
 
 
 
 grid.arrange(P1, P2, P3, P4, ncol=1)
 
}

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).