R/tRackIT_fun.R

In Nature40/tRackIT: Process data recorded with automatic radio-tracking stations

#' calculate triangulations
#'
#' @description calculates triangulations
#'
#'
#' @author Jannis Gottwald
#'
#'
#' @param path_to_directory string, path to directory where files with bearings are stored
#' @param animal list, list generated by initAnimal functiom
#' @param projList list, list generated by initProject functiom
#' @param tw numeric, time window in which bearings are used for triangulation
#' @param be_col character, colname of bearing column
#'
#'
#' @export
#'
#' @examples
#' #projroot<-paste0(getwd(),"/tRackIT_test_data/")
#' # get project file
#' #p <- getProject(projroot)
#' # get meta data of 1 individual
#' #animal <- getAnimal(projroot, animalID = "woodpecker")
#' #tRackIT(animal = animal, projList = p, be_col = "bearings_filtered", tw = 5, path_to_directory = animal$path$bearings_filtered)
#'
tRackIT <- function(animal, projList, tw = 5, be_col = "bearings_filtered", path_to_directory) {

  # define pipe operator

  `%>%` <- dplyr::`%>%`

  # get  hampel filtered data
  fls <- list.files(path_to_directory, full.names = TRUE)

  if (length(fls) == 0) {
    stop("No files available from input directory")
  }

  stations <- projList$stations

  nms_actual <- colnames(stations)
  nms_expected <- c("station", "X", "Y")

  if (!all(nms_expected %in% nms_actual)) {
    idx <- nms_expected %in% nms_actual

    stop(paste0("Required column ", nms_expected[!idx], " in stations not found! "))
  }

  epsg <- projList$epsg

  # read hampel filtered bearings and match station coordinates
  be <- plyr::ldply(fls, function(f) {
    tmp <- data.table::fread(f)

    tmp$y <- stations$Y[stations$station == tmp$station[1]][1]
    tmp$x <- stations$X[stations$station == tmp$station[1]][1]

    return(tmp)
  })


  nms_actual <- colnames(be)

  nms_expected <- c("timestamp", "0", "90", "180", "270", be_col)

  if (!all(nms_expected %in% nms_actual)) {
    idx <- nms_expected %in% nms_actual

    stop(paste0("Required column ", nms_expected[!idx], " in stations not found! "))
  }



  # number of antennas available for bearing calculation
  be$antennas <- rowSums(!is.na(be[, c("0", "90", "180", "270")]))

  # maximum signal strength per bearing
  be$max_dB <- apply(be[, c("0", "90", "180", "270")], 1, max, na.rm = TRUE)

  # cut to timeslots

  be$tc <- cut(be$timestamp, breaks = paste0(tw, " sec"))

  be <- be %>%
    dplyr::group_by(tc) %>%
    dplyr::filter(length(unique(station)) > 1)


  ######### Triangulate
  print("triangulate")
  tri_points <- plyr::ldply(unique(be$tc), .progress = "text", function(ti) {
    
    tmp <- be[be$tc == ti, ]
    # print(head(tmp))

    tmp <- as.data.frame(tmp)
    # print(tmp)

    tryCatch(
      expr = {
        tmp$be <- tmp[, be_col]
        tmp <- as.data.frame(tmp)
        tmp <- tmp[!is.na(tmp$be), ]

        # calculate triangulations with simple intersection method
        tmp2 <- time_match_station(tmp, tmstmp = ti, method = be_col)

        tmp2 <- tri_in(data = tmp2, projList = projList)

        # print(head(tmp2))

        # gather information
        stats <- unique(c(tmp2$name_s1, tmp2$name_s2))
        tmp2$method <- "intersection"

        # calculate triangulations with telemetr package
        sp::coordinates(tmp) <- c("x", "y")
        sp::proj4string(tmp) <- sp::CRS(paste0("+init=epsg:", epsg)) # WGS 84

        tlm <- calc_telemetr_tRackIT(data = tmp, tmstmp = ti)




        # finalize data frame
        df <- as.data.frame(tlm)

        df <- rbind(df, data.frame(ID = df$ID[1], x = tmp2$x, y = tmp2$y, time = ti, sd = NA, kappa = NA, cor = NA, se.x = NA, se.y = NA, ijob = NA, err = NA, npts = NA, method = "intersect"))

        df$number_of_stations <- length(unique(tmp$station))

        df$number_of_potential_biangulations <- nrow(tmp2)

        df$number_of_successfull_biangulations <- nrow(tmp2[!is.na(tmp2$x), ])

        df$min_intersection <- min(tmp2[!is.na(tmp2$x), ]$intersection)

        df$max_intersection <- max(tmp2[!is.na(tmp2$x), ]$intersection)

        df$max_DB <- max(tmp$max_dB)

        df$max_station <- tmp$station[tmp$max_dB == max(tmp$max_dB)][1]



        df$x_max_station <- stations$X[stations$station == df$max_station[1]][1]
        df$y_max_station <- stations$Y[stations$station == df$max_station[1]][1]

        gc()
        # print(df)

        if (is.data.frame(df)) {
          return(df)
        }

        # print(df)
      },
      error = function(e) {
        #   #print(e)
      }
    )
    # print("start gc")

    # print("end gc")
  })
  
  print("calc dist")
  
  print(head(tri_points))
  
  tri_points$dist_to_max_station <- round(raster::pointDistance(cbind(tri_points$x_max_station, tri_points$y_max_station), cbind(tri_points$x, tri_points$y), lonlat = T))
  
  
  
  closest_station <- vector()
  dist_to_closest_station <- vector()
  ant_tst <- stations
  
  print("dist to stations")
  for (n in 1:nrow(tri_points)) {
    
    
    # print(n)
    ant_tst$dist <- NA
    # print(i)
    ant_tst$dist <- round(raster::pointDistance(cbind(tri_points$x[n], tri_points$y[n]), cbind(stations$X, stations$Y), lonlat = T))
    # print(min(ant_tst$dist))
    closest_station[n] <- as.character(ant_tst$station[ant_tst$dist == min(ant_tst$dist)][1])
    
    dist_to_closest_station[n] <- min(ant_tst$dist)
  }
  
  tri_points$closest_station <- closest_station
  tri_points$dist_to_closest_station <- dist_to_closest_station
  
  
  

  data.table::fwrite(tri_points, paste0(animal$path$triangulations, "/triangulations", "_time_window", "_", tw, ".csv"))
  
}
  
 
# dist_to_stats<-function(tri_points){
# 
# 
#  print("calc dist")
# 
#   print(head(tri_points))
# 
#   tri_points$dist_to_max_station <- round(raster::pointDistance(cbind(tri_points$x_max_station, tri_points$y_max_station), cbind(tri_points$x, tri_points$y), lonlat = T))
# 
# 
# 
#   closest_station <- vector()
#   dist_to_closest_station <- vector()
#   ant_tst <- stations
# 
#   print("dist to stations")
#   for (n in 1:nrow(tri_points)) {
# 
# 
#     # print(n)
#     ant_tst$dist <- NA
#     # print(i)
#     ant_tst$dist <- round(raster::pointDistance(cbind(tri_points$x[n], tri_points$y[n]), cbind(stations$X, stations$Y), lonlat = T))
#     # print(min(ant_tst$dist))
#     closest_station[n] <- as.character(ant_tst$station[ant_tst$dist == min(ant_tst$dist)][1])
# 
#     dist_to_closest_station[n] <- min(ant_tst$dist)
#   }
# 
#   tri_points$closest_station <- closest_station
#   tri_points$dist_to_closest_station <- dist_to_closest_station
# 
# 
# 
# 
#   
# }