From 980a4582d3d7ad32fe60251daf92d2db74a463e2 Mon Sep 17 00:00:00 2001 From: Mike O'Brien Date: Thu, 23 May 2024 18:11:28 -0400 Subject: [PATCH 1/8] rm geosphere --- R/vis-interpolate_path.r | 22 +++++++++++++++---- .../detection_range_handout.Rmd | 6 ++--- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/R/vis-interpolate_path.r b/R/vis-interpolate_path.r index 073287ad..1cd8a0ea 100644 --- a/R/vis-interpolate_path.r +++ b/R/vis-interpolate_path.r @@ -337,12 +337,17 @@ interpolate_path <- function(det, trans = NULL, start_time = NULL, on = .(start_dtc >= start, start_dtc <= end)] # calculate great circle distance between coords - #dtc[, gcd := geosphere::distHaversine(as.matrix( + # dtc[, gcd := geosphere::distHaversine(as.matrix( # .SD[1, c("deploy_long", "deploy_lat")]), # as.matrix(.SD[.N, c("deploy_long", "deploy_lat")])), by = i.start] # remove geosphere dependency and use geodist. geosphere relies on sp and raster. - dtc[, gcd := geodist::geodist_vec(x1 = .SD[[1, "deploy_long"]], y1 = .SD[[1, "deploy_lat"]], x2 = .SD[[.N, "deploy_long"]], y2 = .SD[[.N, "deploy_lat"]], measure = "haversine"), by = i.start] + dtc[, gcd := geodist::geodist_vec(x1 = .SD[[1, "deploy_long"]], + y1 = .SD[[1, "deploy_lat"]], + x2 = .SD[[.N, "deploy_long"]], + y2 = .SD[[.N, "deploy_lat"]], + measure = "haversine"), + by = i.start] # calculate least cost (non-linear) distance between points message("Calculating least-cost (non-linear) distances... (step 1 of 3)") @@ -498,8 +503,17 @@ interpolate_path <- function(det, trans = NULL, start_time = NULL, nln_small[, latitude_lead := data.table::shift(nln_latitude, type = "lag", fill = NA), by = i.start] nln_small[, longitude_lead := data.table::shift(nln_longitude, type = "lag", fill = NA), by = i.start] - nln_small[, cumdist := geosphere::distGeo(.SD[, c("nln_longitude", "nln_latitude")], - .SD[,c("longitude_lead", "latitude_lead")]), by = i.start] + # Switch from geosphere to geodist + # nln_small[, cumdist := geosphere::distGeo( + # .SD[, c("nln_longitude", "nln_latitude")], + # .SD[,c("longitude_lead", "latitude_lead")]), by = i.start] + + nln_small[, cumdist := geodist::geodist_vec(x1 = .SD[["nln_longitude"]], + y1 = .SD[["nln_latitude"]], + x2 = .SD[["longitude_lead"]], + y2 = .SD[["latitude_lead"]], + measure = 'geodesic'), + by = .I] nln_small[is.na(cumdist), cumdist := 0] nln_small[, cumdist := cumsum(cumdist), by = i.start] diff --git a/inst/supplemental_docs/detection_range_handout.Rmd b/inst/supplemental_docs/detection_range_handout.Rmd index edb12986..8d192ba6 100644 --- a/inst/supplemental_docs/detection_range_handout.Rmd +++ b/inst/supplemental_docs/detection_range_handout.Rmd @@ -167,7 +167,6 @@ https://github.com/ocean-tracking-network/glatos. # load libraries library(glatos) library(data.table) -library(geosphere) library(ggplot2) # load detection data from glatos package @@ -188,7 +187,7 @@ of columns. Many of the columns in the dtc object are not needed for this analysis. In the next code chunk, we will extract only the necessary columns from the `dtc` object and calculate the geographic distance -between each receiver using the `geosphere::distm` function. The +between each receiver using the `geodist::geodist_vec` function. The internal tag on 4 receivers was used in this study so distances among receivers is the same as distances between each tag and receiver. @@ -202,8 +201,7 @@ dtc <- dtc[, c("detection_timestamp_utc", "transmitter_id", "receiver_sn", rec_dist <- unique(dtc[,c("receiver_sn", "deploy_lat", "deploy_long")]) # calculate distance matrix -dist_matrix <- distm(rec_dist[, c("deploy_long", "deploy_lat")], - rec_dist[, c("deploy_long", "deploy_lat")]) +dist_matrix <- geodist::geodist_vec(rec_dist$deploy_long, rec_dist$deploy_lat) # rename columns and rows of matrix colnames(dist_matrix) <- rec_dist[,get("receiver_sn")] From 28550d60566f1fb4f6813ad53c63b726b2de022b Mon Sep 17 00:00:00 2001 From: Mike O'Brien Date: Fri, 24 May 2024 12:22:21 -0400 Subject: [PATCH 2/8] drop geosphere in favor of hard-coded haversine --- R/util-point_offset.r | 29 ++++++++++++++++++++++++----- man/point_offset.Rd | 3 +-- 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/R/util-point_offset.r b/R/util-point_offset.r index 5c7748fc..7c43db85 100755 --- a/R/util-point_offset.r +++ b/R/util-point_offset.r @@ -1,8 +1,7 @@ #' @title Identify new location based on distance and bearing from another #' #' @description Calculates latitude and longitude for new point that is x meters -#' away at bearing y from a geographic location (Longitude, Latitude). uses -#' "destPoint" function from "geosphere" package and calculations are based on +#' away at bearing y from a geographic location (Longitude, Latitude) using #' great circle distances. #' #' @param lon vector of longitudes (dd) to calculate offset points @@ -32,9 +31,8 @@ point_offset <- function(lon = NA, lat = NA, offsetDist = NA, offsetDir = NA, distUnit = "m") { - # require(geosphere) #for spatialPoints - if (distUnit == "ft") offsetDist <- 0.3048 * offsetDist # conver to m if needed + if (distUnit == "ft") offsetDist <- 0.3048 * offsetDist # convert to m if needed if (!(distUnit) %in% c("ft", "m")) { stop("Input attribute 'dirUnit' must be 'm' (meters) or 'ft' (feet).") } @@ -52,7 +50,28 @@ point_offset <- function(lon = NA, lat = NA, offsetDist = NA, offsetDir = NA, ) bearing <- dirKey$deg[match(offsetDir, dirKey$txt)] + + haversine <- function(lon, lat, bearing, offsetDist){ + lat <- lat * pi/180 + lon <- lon * pi/180 + R <- 6378137 + bearing <- bearing * pi/180 + lat2 <- asin(sin(lat) * cos(offsetDist/R) + cos(lat) * sin(offsetDist/R) * cos(bearing)) + lon2 <- 180/pi * ( + lon + atan2( + sin(bearing) * sin(offsetDist/R) * cos(lat), + cos(offsetDist/R) - sin(lat) * sin(lat2) + ) + ) + + lat2 <- 180/pi * lat2 + + coords <- matrix(c(lon2, lat2), ncol = 2) + colnames(coords) <- c('lon', 'lat') + + return(coords) + } - pos <- geosphere::destPoint(cbind(lon, lat), bearing, offsetDist) + pos <- haversine(lon, lat, bearing, offsetDist) return(pos) } diff --git a/man/point_offset.Rd b/man/point_offset.Rd index de2ce249..494dda33 100644 --- a/man/point_offset.Rd +++ b/man/point_offset.Rd @@ -28,8 +28,7 @@ Options are NA,"N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", } \description{ Calculates latitude and longitude for new point that is x meters -away at bearing y from a geographic location (Longitude, Latitude). uses -"destPoint" function from "geosphere" package and calculations are based on +away at bearing y from a geographic location (Longitude, Latitude) using great circle distances. } \examples{ From dec4ece0dc922abcc2a08fa9eed12eb4e451aa0f Mon Sep 17 00:00:00 2001 From: Mike O'Brien Date: Fri, 24 May 2024 12:22:31 -0400 Subject: [PATCH 3/8] rm geosphere --- DESCRIPTION | 1 - 1 file changed, 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7dc6ebb3..c9b335b1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,7 +13,6 @@ Imports: gdalUtilities, geodist, gdistance, - geosphere, jsonlite, lubridate, magrittr, From f88fcb483024b9befc6176d3332263154e25baeb Mon Sep 17 00:00:00 2001 From: Mike O'Brien Date: Fri, 24 May 2024 12:22:31 -0400 Subject: [PATCH 4/8] rm geosphere. closes #235 --- DESCRIPTION | 1 - 1 file changed, 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7dc6ebb3..c9b335b1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,7 +13,6 @@ Imports: gdalUtilities, geodist, gdistance, - geosphere, jsonlite, lubridate, magrittr, From 64237d6a4bd6d3f2b1e16e28335420a68333b974 Mon Sep 17 00:00:00 2001 From: Mike O'Brien Date: Fri, 24 May 2024 13:23:36 -0400 Subject: [PATCH 5/8] Add initial tests for point_offset --- R/util-point_offset.r | 7 ++- tests/testthat/test-point_offset.R | 69 ++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+), 2 deletions(-) create mode 100644 tests/testthat/test-point_offset.R diff --git a/R/util-point_offset.r b/R/util-point_offset.r index 7c43db85..ee9932bd 100755 --- a/R/util-point_offset.r +++ b/R/util-point_offset.r @@ -66,8 +66,11 @@ point_offset <- function(lon = NA, lat = NA, offsetDist = NA, offsetDir = NA, lat2 <- 180/pi * lat2 - coords <- matrix(c(lon2, lat2), ncol = 2) - colnames(coords) <- c('lon', 'lat') + coords <- matrix( + c(lon2, lat2), + ncol = 2, + dimnames = list(NULL, c('lon', 'lat')) + ) return(coords) } diff --git a/tests/testthat/test-point_offset.R b/tests/testthat/test-point_offset.R new file mode 100644 index 00000000..db02ca93 --- /dev/null +++ b/tests/testthat/test-point_offset.R @@ -0,0 +1,69 @@ +test_that("Returns named matrix", { + expect_type(point_offset(lon = -83.0, lat = 44.0), 'double') + expect_vector(point_offset(lon = -83.0, lat = 44.0)) + + # No inputs + expect_type(point_offset(), 'double') + + expect_identical( + colnames(point_offset(lon = -83.0, lat = 44.0)), + c('lon', 'lat') + ) +}) + + + + +test_that("Converts feet to meters", { + expect_equal( + point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = 'NE', distUnit = 'ft'), + matrix(c(-82.99973, 44.00019), + ncol = 2, + dimnames = list(NULL, c('lon', 'lat'))), + tolerance = 1e-5 + ) + + expect_false( + identical( + point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = 'NE', distUnit = 'm'), + point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = 'NE', distUnit = 'ft') + ) + ) +}) + + + + +test_that("Errors with wrong units", { + expect_error( + point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = 'NE', distUnit = 'km'), + "Input attribute 'dirUnit' must be 'm' \\(meters\\) or 'ft' \\(feet\\)\\." + ) +}) + + + + +test_that("No input returns NA", { + # No distance + expect_true( + all( + is.na( + point_offset(lon = -83.0, lat = 44.0, offsetDir = 'NE') + ) + ) + ) + + # No direction + expect_true( + all( + is.na( + point_offset(lon = -83.0, lat = 44.0, offsetDist = 100) + ) + ) + ) +}) From 0750b5768419892b9573b75e3db56a873508fe6e Mon Sep 17 00:00:00 2001 From: Mike O'Brien Date: Fri, 24 May 2024 14:27:44 -0400 Subject: [PATCH 6/8] Add initial interpolate_path tests --- R/vis-interpolate_path.r | 5 +- tests/testthat/_snaps/interpolate_path.md | 136 +++++++++++++ tests/testthat/test-interpolate_path.R | 227 ++++++++++++++++++++++ 3 files changed, 367 insertions(+), 1 deletion(-) create mode 100644 tests/testthat/_snaps/interpolate_path.md create mode 100644 tests/testthat/test-interpolate_path.R diff --git a/R/vis-interpolate_path.r b/R/vis-interpolate_path.r index 1cd8a0ea..a77ad91a 100644 --- a/R/vis-interpolate_path.r +++ b/R/vis-interpolate_path.r @@ -199,7 +199,10 @@ interpolate_path <- function(det, trans = NULL, start_time = NULL, # stop if out_class is not NULL, data.table, or tibble if(!is.null(out_class)){ - if( !(out_class %in% c("data.table", "tibble"))) {stop('out_class is not a "data.table" or "tibble"')}} + if( !(out_class %in% c("data.table", "tibble"))) { + stop('out_class is not a "data.table" or "tibble"') + } + } # check to see that trans is a transition layer or transition stack if(!is.null(trans) & diff --git a/tests/testthat/_snaps/interpolate_path.md b/tests/testthat/_snaps/interpolate_path.md new file mode 100644 index 00000000..cb856336 --- /dev/null +++ b/tests/testthat/_snaps/interpolate_path.md @@ -0,0 +1,136 @@ +# linear interpolation works + + Code + linear_interp + Output + animal_id bin_timestamp latitude longitude record_type + 1 1 2000-01-01 44.00000 -87.00000 detection + 2 1 2000-01-02 44.01613 -86.85484 interpolated + 3 1 2000-01-03 44.03226 -86.70968 interpolated + 4 1 2000-01-04 44.04839 -86.56452 interpolated + 5 1 2000-01-05 44.06452 -86.41935 interpolated + 6 1 2000-01-06 44.08065 -86.27419 interpolated + 7 1 2000-01-07 44.09677 -86.12903 interpolated + 8 1 2000-01-08 44.11290 -85.98387 interpolated + 9 1 2000-01-09 44.12903 -85.83871 interpolated + 10 1 2000-01-10 44.14516 -85.69355 interpolated + 11 1 2000-01-11 44.16129 -85.54839 interpolated + 12 1 2000-01-12 44.17742 -85.40323 interpolated + 13 1 2000-01-13 44.19355 -85.25806 interpolated + 14 1 2000-01-14 44.20968 -85.11290 interpolated + 15 1 2000-01-15 44.22581 -84.96774 interpolated + 16 1 2000-01-16 44.24194 -84.82258 interpolated + 17 1 2000-01-17 44.25806 -84.67742 interpolated + 18 1 2000-01-18 44.27419 -84.53226 interpolated + 19 1 2000-01-19 44.29032 -84.38710 interpolated + 20 1 2000-01-20 44.30645 -84.24194 interpolated + 21 1 2000-01-21 44.32258 -84.09677 interpolated + 22 1 2000-01-22 44.33871 -83.95161 interpolated + 23 1 2000-01-23 44.35484 -83.80645 interpolated + 24 1 2000-01-24 44.37097 -83.66129 interpolated + 25 1 2000-01-25 44.38710 -83.51613 interpolated + 26 1 2000-01-26 44.40323 -83.37097 interpolated + 27 1 2000-01-27 44.41935 -83.22581 interpolated + 28 1 2000-01-28 44.43548 -83.08065 interpolated + 29 1 2000-01-29 44.45161 -82.93548 interpolated + 30 1 2000-01-30 44.46774 -82.79032 interpolated + 31 1 2000-01-31 44.48387 -82.64516 interpolated + 32 1 2000-02-01 44.50000 -82.50000 detection + 33 1 2000-02-02 44.46552 -82.34483 interpolated + 34 1 2000-02-03 44.43103 -82.18966 interpolated + 35 1 2000-02-04 44.39655 -82.03448 interpolated + 36 1 2000-02-05 44.36207 -81.87931 interpolated + 37 1 2000-02-06 44.32759 -81.72414 interpolated + 38 1 2000-02-07 44.29310 -81.56897 interpolated + 39 1 2000-02-08 44.25862 -81.41379 interpolated + 40 1 2000-02-09 44.22414 -81.25862 interpolated + 41 1 2000-02-10 44.18966 -81.10345 interpolated + 42 1 2000-02-11 44.15517 -80.94828 interpolated + 43 1 2000-02-12 44.12069 -80.79310 interpolated + 44 1 2000-02-13 44.08621 -80.63793 interpolated + 45 1 2000-02-14 44.05172 -80.48276 interpolated + 46 1 2000-02-15 44.01724 -80.32759 interpolated + 47 1 2000-02-16 43.98276 -80.17241 interpolated + 48 1 2000-02-17 43.94828 -80.01724 interpolated + 49 1 2000-02-18 43.91379 -79.86207 interpolated + 50 1 2000-02-19 43.87931 -79.70690 interpolated + 51 1 2000-02-20 43.84483 -79.55172 interpolated + 52 1 2000-02-21 43.81034 -79.39655 interpolated + 53 1 2000-02-22 43.77586 -79.24138 interpolated + 54 1 2000-02-23 43.74138 -79.08621 interpolated + 55 1 2000-02-24 43.70690 -78.93103 interpolated + 56 1 2000-02-25 43.67241 -78.77586 interpolated + 57 1 2000-02-26 43.63793 -78.62069 interpolated + 58 1 2000-02-27 43.60345 -78.46552 interpolated + 59 1 2000-02-28 43.56897 -78.31034 interpolated + 60 1 2000-02-29 43.53448 -78.15517 interpolated + 61 1 2000-03-01 43.50000 -78.00000 detection + +# Non-linear interpolation works + + Code + nonlinear_interp + Output + animal_id bin_timestamp latitude longitude record_type + 1 1 2000-01-01 44.00000 -87.00000 detection + 2 1 2000-01-02 44.14456 -87.01280 interpolated + 3 1 2000-01-03 44.28736 -86.99353 interpolated + 4 1 2000-01-04 44.39165 -86.84604 interpolated + 5 1 2000-01-05 44.49603 -86.69840 interpolated + 6 1 2000-01-06 44.60050 -86.55064 interpolated + 7 1 2000-01-07 44.70507 -86.40275 interpolated + 8 1 2000-01-08 44.80973 -86.25472 interpolated + 9 1 2000-01-09 44.91448 -86.10656 interpolated + 10 1 2000-01-10 45.01933 -85.95826 interpolated + 11 1 2000-01-11 45.12427 -85.80983 interpolated + 12 1 2000-01-12 45.22932 -85.66127 interpolated + 13 1 2000-01-13 45.33445 -85.51257 interpolated + 14 1 2000-01-14 45.43968 -85.36373 interpolated + 15 1 2000-01-15 45.54501 -85.21476 interpolated + 16 1 2000-01-16 45.65044 -85.06565 interpolated + 17 1 2000-01-17 45.74961 -84.91262 interpolated + 18 1 2000-01-18 45.75648 -84.70329 interpolated + 19 1 2000-01-19 45.72012 -84.49758 interpolated + 20 1 2000-01-20 45.65750 -84.30603 interpolated + 21 1 2000-01-21 45.59065 -84.11694 interpolated + 22 1 2000-01-22 45.52387 -83.92803 interpolated + 23 1 2000-01-23 45.45715 -83.73930 interpolated + 24 1 2000-01-24 45.37617 -83.56547 interpolated + 25 1 2000-01-25 45.27099 -83.41672 interpolated + 26 1 2000-01-26 45.16591 -83.26810 interpolated + 27 1 2000-01-27 45.06093 -83.11961 interpolated + 28 1 2000-01-28 44.95605 -82.97127 interpolated + 29 1 2000-01-29 44.85125 -82.82305 interpolated + 30 1 2000-01-30 44.74656 -82.67497 interpolated + 31 1 2000-01-31 44.64195 -82.52702 interpolated + 32 1 2000-02-01 44.50000 -82.50000 detection + 33 1 2000-02-02 44.27598 -82.41394 interpolated + 34 1 2000-02-03 44.02721 -82.41394 interpolated + 35 1 2000-02-04 43.77843 -82.41394 interpolated + 36 1 2000-02-05 43.52964 -82.41394 interpolated + 37 1 2000-02-06 43.28084 -82.41394 interpolated + 38 1 2000-02-07 43.03202 -82.41394 interpolated + 39 1 2000-02-08 42.80359 -82.47917 interpolated + 40 1 2000-02-09 42.60977 -82.59843 interpolated + 41 1 2000-02-10 42.43783 -82.84162 interpolated + 42 1 2000-02-11 42.27309 -83.09093 interpolated + 43 1 2000-02-12 42.03710 -83.06659 interpolated + 44 1 2000-02-13 41.99075 -82.75061 interpolated + 45 1 2000-02-14 42.02156 -82.42462 interpolated + 46 1 2000-02-15 42.12825 -82.12284 interpolated + 47 1 2000-02-16 42.23508 -81.82064 interpolated + 48 1 2000-02-17 42.34205 -81.51804 interpolated + 49 1 2000-02-18 42.43982 -81.21055 interpolated + 50 1 2000-02-19 42.49761 -80.88362 interpolated + 51 1 2000-02-20 42.55545 -80.55641 interpolated + 52 1 2000-02-21 42.64289 -80.24305 interpolated + 53 1 2000-02-22 42.71383 -79.92129 interpolated + 54 1 2000-02-23 42.77186 -79.59301 interpolated + 55 1 2000-02-24 42.82994 -79.26443 interpolated + 56 1 2000-02-25 42.88807 -78.93556 interpolated + 57 1 2000-02-26 43.07781 -79.02641 interpolated + 58 1 2000-02-27 43.28679 -78.97607 interpolated + 59 1 2000-02-28 43.38565 -78.66456 interpolated + 60 1 2000-02-29 43.44428 -78.33286 interpolated + 61 1 2000-03-01 43.50000 -78.00000 detection + diff --git a/tests/testthat/test-interpolate_path.R b/tests/testthat/test-interpolate_path.R new file mode 100644 index 00000000..a144bf8d --- /dev/null +++ b/tests/testthat/test-interpolate_path.R @@ -0,0 +1,227 @@ +# Set up data from interpolate_path example 1 +pos <- data.frame( + animal_id=1, + deploy_long=c(-87,-82.5, -78), + deploy_lat=c(44, 44.5, 43.5), + detection_timestamp_utc=as.POSIXct( + c("2000-01-01 00:00", "2000-02-01 00:00", "2000-03-01 00:00"), + tz = "UTC" + ) +) + + + + +test_that("linear interpolation works", { + linear_interp <- interpolate_path(pos) + expect_snapshot( + linear_interp + ) + + expect_s3_class( + linear_interp, + 'data.frame' + ) + + expect_named( + linear_interp, + c("animal_id", "bin_timestamp", "latitude", "longitude", "record_type") + ) + + # Retains actual detections + expect_equal( + nrow( + linear_interp[linear_interp$record_type == 'detection',] + ), + nrow(pos) + ) +}) + + + + +test_that("Non-linear interpolation works", { + linear_interp <- interpolate_path(pos) + expect_message( + nonlinear_interp <- interpolate_path(pos, trans = greatLakesTrLayer), + 'Calculating least-cost \\(non-linear\\) distances\\.\\.\\. \\(step 1 of 3\\)' + ) |> + expect_output('|=*| 100%') |> + expect_message('Starting non-linear interpolation\\.\\.\\. \\(step 3 of 3\\)') |> + expect_message('Finalizing results\\.') + + expect_false( + identical( + linear_interp, + nonlinear_interp + ) + ) + + expect_snapshot( + nonlinear_interp + ) + + expect_s3_class( + nonlinear_interp, + 'data.frame' + ) + + expect_named( + nonlinear_interp, + c("animal_id", "bin_timestamp", "latitude", "longitude", "record_type") + ) + + # Retains actual detections + expect_equal( + nrow( + nonlinear_interp[nonlinear_interp$record_type == 'detection',] + ), + nrow(pos) + ) +}) + + + +test_that("Forced linear interpolation works", { + linear_interp <- interpolate_path(pos) + nonlinear_interp <- suppressMessages( + interpolate_path(pos, trans = greatLakesTrLayer, + show_progress = FALSE) + ) + + + forced_linear_interp <- interpolate_path(pos, trans = greatLakesTrLayer, + lnl_thresh = 0) |> + expect_no_message() + + + expect_false( + identical( + forced_linear_interp, + nonlinear_interp + ) + ) + + expect_identical( + linear_interp, + forced_linear_interp + ) + + expect_s3_class( + forced_linear_interp, + 'data.frame' + ) + + expect_named( + forced_linear_interp, + c("animal_id", "bin_timestamp", "latitude", "longitude", "record_type") + ) + + # Retains actual detections + expect_equal( + nrow( + forced_linear_interp[forced_linear_interp$record_type == 'detection',] + ), + nrow(pos) + ) +}) + + + +test_that("Progress bar can be suppressed", { + expect_output( + suppressMessages( + interpolate_path(pos, trans = greatLakesTrLayer, + show_progress = FALSE) + ), + NA + ) +}) + + + + +test_that("Checks output class", { + # Error if wrong class + expect_error( + interpolate_path(pos, out_class = 'list'), + "out_class is not a \"data\\.table\" or \"tibble\"" + ) + + # data.frame if NULL + expect_s3_class( + interpolate_path(pos, out_class = NULL), + 'data.frame', + exact = TRUE + ) + + # data.table + expect_s3_class( + interpolate_path(pos, out_class = 'data.table'), + c('data.table', 'data.frame'), + exact = TRUE + ) + + # tibble + expect_s3_class( + interpolate_path(pos, out_class = 'tibble'), + c('tbl_df', 'tbl', 'data.frame'), + exact = TRUE + ) +}) + + + +test_that("Checks trans is a transition layer or transition stack", { + suppressMessages( + trans <- make_transition3(great_lakes_polygon) + ) + + expect_error( + interpolate_path(pos, trans = trans), + "Supplied object for 'trans' argument is not class TransitionLayer or TransitionStack\\." + ) +}) + + + +test_that("Checks that start_time was successfully coverted", { + expect_error( + interpolate_path(pos, start_time = NA), + "start_time cannot be coerced to 'POSIXct' or 'POSIXt' class" + ) +}) + + + +test_that("Checks that start_time < largest timestamp in dataset", { + expect_error( + interpolate_path(pos, start_time = max(pos$detection_timestamp_utc) + 1), + "start_time is larger than last detection\\. No data to interpolate\\!" + ) +}) + + + +test_that("Checks that fish have multiple observations", { + expect_error( + interpolate_path(pos[1,]), + "must have two observations to interpolate" + ) +}) + + + +test_that("Errors and displays offending receivers if any receivers are on land", { + pos2 <- pos + pos2$deploy_long[1] <- -90 + + expect_error( + suppressMessages( + interpolate_path(pos2, trans = greatLakesTrLayer, show_progress = FALSE) + ), + "Some coordinates are on land or beyond extent. + Interpolation impossible! Check receiver locations or extents of transition + layer:\n" + ) +}) From c64e8c9a905f3c7b4bf8d898cff81231172f071a Mon Sep 17 00:00:00 2001 From: mhpob Date: Fri, 24 May 2024 18:30:35 +0000 Subject: [PATCH 7/8] Style code (GHA) --- R/util-point_offset.r | 29 ++++--- tests/testthat/test-interpolate_path.R | 104 +++++++++++++------------ tests/testthat/test-point_offset.R | 45 ++++++----- 3 files changed, 96 insertions(+), 82 deletions(-) diff --git a/R/util-point_offset.r b/R/util-point_offset.r index ee9932bd..5a255777 100755 --- a/R/util-point_offset.r +++ b/R/util-point_offset.r @@ -31,7 +31,6 @@ point_offset <- function(lon = NA, lat = NA, offsetDist = NA, offsetDir = NA, distUnit = "m") { - if (distUnit == "ft") offsetDist <- 0.3048 * offsetDist # convert to m if needed if (!(distUnit) %in% c("ft", "m")) { stop("Input attribute 'dirUnit' must be 'm' (meters) or 'ft' (feet).") @@ -50,28 +49,28 @@ point_offset <- function(lon = NA, lat = NA, offsetDist = NA, offsetDir = NA, ) bearing <- dirKey$deg[match(offsetDir, dirKey$txt)] - - haversine <- function(lon, lat, bearing, offsetDist){ - lat <- lat * pi/180 - lon <- lon * pi/180 + + haversine <- function(lon, lat, bearing, offsetDist) { + lat <- lat * pi / 180 + lon <- lon * pi / 180 R <- 6378137 - bearing <- bearing * pi/180 - lat2 <- asin(sin(lat) * cos(offsetDist/R) + cos(lat) * sin(offsetDist/R) * cos(bearing)) - lon2 <- 180/pi * ( + bearing <- bearing * pi / 180 + lat2 <- asin(sin(lat) * cos(offsetDist / R) + cos(lat) * sin(offsetDist / R) * cos(bearing)) + lon2 <- 180 / pi * ( lon + atan2( - sin(bearing) * sin(offsetDist/R) * cos(lat), - cos(offsetDist/R) - sin(lat) * sin(lat2) + sin(bearing) * sin(offsetDist / R) * cos(lat), + cos(offsetDist / R) - sin(lat) * sin(lat2) ) ) - - lat2 <- 180/pi * lat2 - + + lat2 <- 180 / pi * lat2 + coords <- matrix( c(lon2, lat2), ncol = 2, - dimnames = list(NULL, c('lon', 'lat')) + dimnames = list(NULL, c("lon", "lat")) ) - + return(coords) } diff --git a/tests/testthat/test-interpolate_path.R b/tests/testthat/test-interpolate_path.R index a144bf8d..50ce3dbd 100644 --- a/tests/testthat/test-interpolate_path.R +++ b/tests/testthat/test-interpolate_path.R @@ -1,9 +1,9 @@ # Set up data from interpolate_path example 1 pos <- data.frame( - animal_id=1, - deploy_long=c(-87,-82.5, -78), - deploy_lat=c(44, 44.5, 43.5), - detection_timestamp_utc=as.POSIXct( + animal_id = 1, + deploy_long = c(-87, -82.5, -78), + deploy_lat = c(44, 44.5, 43.5), + detection_timestamp_utc = as.POSIXct( c("2000-01-01 00:00", "2000-02-01 00:00", "2000-03-01 00:00"), tz = "UTC" ) @@ -17,21 +17,21 @@ test_that("linear interpolation works", { expect_snapshot( linear_interp ) - + expect_s3_class( linear_interp, - 'data.frame' + "data.frame" ) - + expect_named( linear_interp, c("animal_id", "bin_timestamp", "latitude", "longitude", "record_type") ) - + # Retains actual detections expect_equal( nrow( - linear_interp[linear_interp$record_type == 'detection',] + linear_interp[linear_interp$record_type == "detection", ] ), nrow(pos) ) @@ -44,37 +44,37 @@ test_that("Non-linear interpolation works", { linear_interp <- interpolate_path(pos) expect_message( nonlinear_interp <- interpolate_path(pos, trans = greatLakesTrLayer), - 'Calculating least-cost \\(non-linear\\) distances\\.\\.\\. \\(step 1 of 3\\)' - ) |> - expect_output('|=*| 100%') |> - expect_message('Starting non-linear interpolation\\.\\.\\. \\(step 3 of 3\\)') |> - expect_message('Finalizing results\\.') - + "Calculating least-cost \\(non-linear\\) distances\\.\\.\\. \\(step 1 of 3\\)" + ) |> + expect_output("|=*| 100%") |> + expect_message("Starting non-linear interpolation\\.\\.\\. \\(step 3 of 3\\)") |> + expect_message("Finalizing results\\.") + expect_false( identical( linear_interp, nonlinear_interp ) ) - + expect_snapshot( nonlinear_interp ) - + expect_s3_class( nonlinear_interp, - 'data.frame' + "data.frame" ) - + expect_named( nonlinear_interp, c("animal_id", "bin_timestamp", "latitude", "longitude", "record_type") ) - + # Retains actual detections expect_equal( nrow( - nonlinear_interp[nonlinear_interp$record_type == 'detection',] + nonlinear_interp[nonlinear_interp$record_type == "detection", ] ), nrow(pos) ) @@ -85,42 +85,46 @@ test_that("Non-linear interpolation works", { test_that("Forced linear interpolation works", { linear_interp <- interpolate_path(pos) nonlinear_interp <- suppressMessages( - interpolate_path(pos, trans = greatLakesTrLayer, - show_progress = FALSE) + interpolate_path(pos, + trans = greatLakesTrLayer, + show_progress = FALSE + ) ) - - - forced_linear_interp <- interpolate_path(pos, trans = greatLakesTrLayer, - lnl_thresh = 0) |> + + + forced_linear_interp <- interpolate_path(pos, + trans = greatLakesTrLayer, + lnl_thresh = 0 + ) |> expect_no_message() - - + + expect_false( identical( forced_linear_interp, nonlinear_interp ) ) - + expect_identical( linear_interp, forced_linear_interp ) - + expect_s3_class( forced_linear_interp, - 'data.frame' + "data.frame" ) - + expect_named( forced_linear_interp, c("animal_id", "bin_timestamp", "latitude", "longitude", "record_type") ) - + # Retains actual detections expect_equal( nrow( - forced_linear_interp[forced_linear_interp$record_type == 'detection',] + forced_linear_interp[forced_linear_interp$record_type == "detection", ] ), nrow(pos) ) @@ -131,8 +135,10 @@ test_that("Forced linear interpolation works", { test_that("Progress bar can be suppressed", { expect_output( suppressMessages( - interpolate_path(pos, trans = greatLakesTrLayer, - show_progress = FALSE) + interpolate_path(pos, + trans = greatLakesTrLayer, + show_progress = FALSE + ) ), NA ) @@ -144,28 +150,28 @@ test_that("Progress bar can be suppressed", { test_that("Checks output class", { # Error if wrong class expect_error( - interpolate_path(pos, out_class = 'list'), + interpolate_path(pos, out_class = "list"), "out_class is not a \"data\\.table\" or \"tibble\"" ) - + # data.frame if NULL expect_s3_class( interpolate_path(pos, out_class = NULL), - 'data.frame', + "data.frame", exact = TRUE ) - + # data.table expect_s3_class( - interpolate_path(pos, out_class = 'data.table'), - c('data.table', 'data.frame'), + interpolate_path(pos, out_class = "data.table"), + c("data.table", "data.frame"), exact = TRUE ) - + # tibble expect_s3_class( - interpolate_path(pos, out_class = 'tibble'), - c('tbl_df', 'tbl', 'data.frame'), + interpolate_path(pos, out_class = "tibble"), + c("tbl_df", "tbl", "data.frame"), exact = TRUE ) }) @@ -176,7 +182,7 @@ test_that("Checks trans is a transition layer or transition stack", { suppressMessages( trans <- make_transition3(great_lakes_polygon) ) - + expect_error( interpolate_path(pos, trans = trans), "Supplied object for 'trans' argument is not class TransitionLayer or TransitionStack\\." @@ -205,7 +211,7 @@ test_that("Checks that start_time < largest timestamp in dataset", { test_that("Checks that fish have multiple observations", { expect_error( - interpolate_path(pos[1,]), + interpolate_path(pos[1, ]), "must have two observations to interpolate" ) }) @@ -215,7 +221,7 @@ test_that("Checks that fish have multiple observations", { test_that("Errors and displays offending receivers if any receivers are on land", { pos2 <- pos pos2$deploy_long[1] <- -90 - + expect_error( suppressMessages( interpolate_path(pos2, trans = greatLakesTrLayer, show_progress = FALSE) diff --git a/tests/testthat/test-point_offset.R b/tests/testthat/test-point_offset.R index db02ca93..83de1a93 100644 --- a/tests/testthat/test-point_offset.R +++ b/tests/testthat/test-point_offset.R @@ -1,13 +1,13 @@ test_that("Returns named matrix", { - expect_type(point_offset(lon = -83.0, lat = 44.0), 'double') + expect_type(point_offset(lon = -83.0, lat = 44.0), "double") expect_vector(point_offset(lon = -83.0, lat = 44.0)) - + # No inputs - expect_type(point_offset(), 'double') - + expect_type(point_offset(), "double") + expect_identical( colnames(point_offset(lon = -83.0, lat = 44.0)), - c('lon', 'lat') + c("lon", "lat") ) }) @@ -16,20 +16,27 @@ test_that("Returns named matrix", { test_that("Converts feet to meters", { expect_equal( - point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, - offsetDir = 'NE', distUnit = 'ft'), + point_offset( + lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = "NE", distUnit = "ft" + ), matrix(c(-82.99973, 44.00019), - ncol = 2, - dimnames = list(NULL, c('lon', 'lat'))), + ncol = 2, + dimnames = list(NULL, c("lon", "lat")) + ), tolerance = 1e-5 ) - + expect_false( identical( - point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, - offsetDir = 'NE', distUnit = 'm'), - point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, - offsetDir = 'NE', distUnit = 'ft') + point_offset( + lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = "NE", distUnit = "m" + ), + point_offset( + lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = "NE", distUnit = "ft" + ) ) ) }) @@ -39,8 +46,10 @@ test_that("Converts feet to meters", { test_that("Errors with wrong units", { expect_error( - point_offset(lon = -83.0, lat = 44.0, offsetDist = 100, - offsetDir = 'NE', distUnit = 'km'), + point_offset( + lon = -83.0, lat = 44.0, offsetDist = 100, + offsetDir = "NE", distUnit = "km" + ), "Input attribute 'dirUnit' must be 'm' \\(meters\\) or 'ft' \\(feet\\)\\." ) }) @@ -53,11 +62,11 @@ test_that("No input returns NA", { expect_true( all( is.na( - point_offset(lon = -83.0, lat = 44.0, offsetDir = 'NE') + point_offset(lon = -83.0, lat = 44.0, offsetDir = "NE") ) ) ) - + # No direction expect_true( all( From d4ca28be989881aa5c6270b20d61b35d056f8cad Mon Sep 17 00:00:00 2001 From: chrisholbrook Date: Fri, 16 Aug 2024 20:51:20 +0000 Subject: [PATCH 8/8] Style code (GHA) --- R/vis-interpolate_path.r | 1406 +++++++++++++++++++------------------- 1 file changed, 706 insertions(+), 700 deletions(-) diff --git a/R/vis-interpolate_path.r b/R/vis-interpolate_path.r index 60ba3731..67e71247 100644 --- a/R/vis-interpolate_path.r +++ b/R/vis-interpolate_path.r @@ -1,700 +1,706 @@ -#' Interpolate new positions within a spatiotemporal path data -#' -#' Interpolate new positions within a spatiotemporal path data set -#' (e.g., detections of tagged fish) at regularly-spaced time intervals -#' using linear or non-linear interpolation. -#' -#' @param det An object of class `glatos_detections` or data frame -#' containing spatiotemporal data with at least 4 columns containing -#' 'animal_id', 'detection_timestamp_utc', 'deploy_long', and -#' 'deploy_lat' columns. -#' -#' @param start_time specify the first time bin for interpolated data. -#' If not supplied, default is first timestamp in the input data -#' set. Must be a character string that can be coerced to -#' 'POSIXct' or an object of class 'POSIXct'. If character string -#' is supplied, timezone is automatically set to UTC. -#' -#' @param out_class Return results as a data.table or tibble. Default -#' returns results as data.frame. Accepts `data.table` or `tibble`. -#' -#' @param int_time_stamp The time step size (in seconds) of interpolated -#' positions. Default is 86400 (one day). -#' -#' @param trans An optional transition matrix with the "cost" of -#' moving across each cell within the map extent. Must be of class -#' `TransitionLayer`. A transition layer may be -#' created from a polygon shapefile using [make_transition]. -#' -#' @param lnl_thresh A numeric threshold for determining if linear or -#' non-linear interpolation shortest path will be used. -#' -#' @param show_progress Logical. Progress bar and status messages will be -#' shown if TRUE (default) and not shown if FALSE. -#' -#' @details Non-linear interpolation uses the `gdistance` package -#' to find the shortest pathway between two locations (i.e., -#' receivers) that avoid 'impossible' movements (e.g., over land for -#' fish). The shortest non-linear path between two locations is -#' calculated using a transition matrix layer that represents the -#' 'cost' of an animal moving between adjacent grid cells. The -#' transition matrix layer (see [gdistance]) is created from -#' a polygon shapefile using [make_transition] or from a -#' `RasterLayer` object using [transition][gdistance::transition]. In -#' `make_transition`, each cell in the output transition layer -#' is coded as water (1) or land (0) to represent possible (1) and -#' impossible (0) movement paths. -#' -#' @details Linear interpolation is used for all points when -#' `trans` is not supplied. When `trans` is supplied, -#' then interpolation method is determined for each pair of -#' sequential observed detections. For example, linear interpolation -#' will be used if the two geographical positions are exactly the -#' same and when the ratio (linear distance:non-linear distance) -#' between two positions is less than `lnl_thresh`. Non-linear -#' interpolation will be used when ratio is greater than -#' `lnl_thresh`. When the ratio of linear distance to -#' non-linear distance is greater than `lnl_thresh`, then the -#' distance of the non-linear path needed to avoid land is greater -#' than the linear path that crosses land. `lnl_thresh` can be -#' used to control whether non-linear or linear interpolation is -#' used for all points. For example, non-linear interpolation will -#' be used for all points when `lnl_thresh` > 1 and linear -#' interpolation will be used for all points when `lnl_thresh` -#' = 0. -#' -#' @details All linear interpolation is done by code{stats::approx} with -#' argument `ties = "ordered"` controlling how tied `x` values -#' are handled. See [approxfun()]. -#' -#' @return A dataframe with animal_id, bin_timestamp, -#' latitude, longitude, and record_type. -#' -#' -#' @author Todd Hayden, Tom Binder, Chris Holbrook -#' -#' @examples -#' -#' #-------------------------------------------------- -#' # EXAMPLE #1 - simple interpolate among lakes -#' -#' # get polygon of the Great Lakes -#' data(great_lakes_polygon) # glatos example data -#' plot(sf::st_geometry(great_lakes_polygon), xlim = c(-92, -76)) -#' -#' # make sample detections data frame -#' pos <- data.frame( -#' animal_id = 1, -#' deploy_long = c(-87, -82.5, -78), -#' deploy_lat = c(44, 44.5, 43.5), -#' detection_timestamp_utc = as.POSIXct(c( -#' "2000-01-01 00:00", -#' "2000-02-01 00:00", "2000-03-01 00:00" -#' ), tz = "UTC") -#' ) -#' -#' # add to plot -#' points(deploy_lat ~ deploy_long, data = pos, pch = 20, cex = 2, col = "red") -#' -#' # interpolate path using linear method -#' path1 <- interpolate_path(pos) -#' nrow(path1) # now 61 points -#' sum(path1$record_type == "interpolated") # 58 interpolated points -#' -#' # add linear path to plot -#' points(latitude ~ longitude, data = path1, pch = 20, cex = 0.8, col = "blue") -#' -#' # load a transition matrix of Great Lakes -#' # NOTE: This is a LOW RESOLUTION TransitionLayer suitable only for -#' # coarse/large scale interpolation only. Most realistic uses -#' # will need to create a TransitionLayer; see ?make_transition. -#' data(greatLakesTrLayer) # glatos example data; a TransitionLayer -#' -#' # interpolate path using non-linear method (requires 'trans') -#' path2 <- interpolate_path(pos, trans = greatLakesTrLayer) -#' -#' # add non-linear path to plot -#' points(latitude ~ longitude, -#' data = path2, pch = 20, cex = 1, -#' col = "green" -#' ) -#' -#' # can also force linear-interpolation with lnlThresh = 0 -#' path3 <- interpolate_path(pos, trans = greatLakesTrLayer, lnl_thresh = 0) -#' -#' # add new linear path to plot -#' points(latitude ~ longitude, -#' data = path3, pch = 20, cex = 1, -#' col = "magenta" -#' ) -#' -#' #-------------------------------------------------- -#' # EXAMPLE #2 - walleye in western Lake Erie -#' \dontrun{ -#' -#' -#' # get example walleye detection data -#' det_file <- system.file("extdata", "walleye_detections.csv", -#' package = "glatos" -#' ) -#' det <- read_glatos_detections(det_file) -#' -#' # take a look -#' head(det) -#' -#' # extract one fish and subset date -#' det <- det[det$animal_id == 22 & -#' det$detection_timestamp_utc > as.POSIXct("2012-04-08") & -#' det$detection_timestamp_utc < as.POSIXct("2013-04-15"), ] -#' -#' # get polygon of the Great Lakes -#' data(great_lakes_polygon) # glatos example data; an sf object -#' -#' # convert polygon to terra::spatVector -#' great_lakes_polygon <- terra::vect(great_lakes_polygon) -#' -#' # crop polygon to western Lake Erie -#' maumee <- terra::crop(great_lakes_polygon, -#' y = terra::ext(-83.7, -82.5, 41.3, 42.4) -#' ) -#' -#' -#' plot(maumee, col = "grey") -#' points(deploy_lat ~ deploy_long, -#' data = det, pch = 20, col = "red", -#' xlim = c(-83.7, -80) -#' ) -#' -#' -#' # make transition layer object -#' tran <- make_transition3(sf::st_as_sf(maumee), res = c(0.1, 0.1)) -#' -#' # plot to check output -#' plot(tran$rast, xlim = c(-83.7, -82.0), ylim = c(41.3, 42.7)) -#' plot(maumee, add = TRUE) -#' -#' # not high enough resolution- bump up resolution, will take some time -#' tran1 <- make_transition3(sf::st_as_sf(maumee), res = c(0.001, 0.001)) -#' -#' # plot to check resolution- much better -#' plot(tran1$rast, xlim = c(-83.7, -82.0), ylim = c(41.3, 42.7)) -#' plot(maumee, add = TRUE) -#' -#' -#' # add fish detections to make sure they are "on the map" -#' # plot unique values only for simplicity -#' foo <- unique(det[, c("deploy_lat", "deploy_long")]) -#' points(foo$deploy_long, foo$deploy_lat, pch = 20, col = "red") -#' -#' # call with "transition matrix" (non-linear interpolation), other options -#' # note that it is quite a bit slower than linear interpolation -#' pos2 <- interpolate_path(det, -#' trans = tran1$transition, -#' out_class = "data.table" -#' ) -#' -#' plot(maumee, col = "grey") -#' points(latitude ~ longitude, data = pos2, pch = 20, col = "red", cex = 0.5) -#' } -#' -#' @export - - - -interpolate_path <- function(det, trans = NULL, start_time = NULL, - int_time_stamp = 86400, lnl_thresh = 0.9, - out_class = NULL, show_progress = TRUE) { - # stop if out_class is not NULL, data.table, or tibble - if (!is.null(out_class)) { - if (!(out_class %in% c("data.table", "tibble"))) { - stop('out_class is not a "data.table" or "tibble"') - } - } - - # check to see that trans is a transition layer or transition stack - if (!is.null(trans) & - inherits(trans, c("TransitionLayer", "TransitionStack")) == FALSE) { - stop( - paste0( - "Supplied object for 'trans' argument is not class ", - "TransitionLayer or TransitionStack." - ), - call. = FALSE - ) - } - - # check start_time - if (is.null(start_time)) { - start_time <- min(det$detection_timestamp_utc) - } - if (is.na(start_time) & length(start_time) > 0) { - stop("start_time cannot be coerced to 'POSIXct' or 'POSIXt' class") - } - - if (is.character(start_time)) { - start_time <- as.POSIXct(start_time, tz = "UTC") - } - - # make sure start_time < largest timestamp in dataset - if (start_time > max(det$detection_timestamp_utc)) { - stop("start_time is larger than last detection. No data to interpolate!", call. = FALSE) - } - - # make copy of detections for function - dtc <- data.table::as.data.table(det) - - # subset only columns for function and rows >= start_time: - dtc <- dtc[detection_timestamp_utc >= start_time, c( - "animal_id", - "detection_timestamp_utc", - "deploy_lat", - "deploy_long" - )] - - dtc[, record_type := "detection"] - - # count number of rows- single observations are not interpolated - dtc[, num_rows := nrow(.SD), by = animal_id] - - # Sort detections by transmitter id and then by detection timestamp - data.table::setkey(dtc, animal_id, detection_timestamp_utc) - - # save original dataset to combine with interpolated data in the end - det <- data.table::copy(dtc) - data.table::setnames(det, c( - "animal_id", "bin_stamp", "i_lat", "i_lon", - "record_type", "num_rows" - )) - - # remove any fish with only one detection - dtc <- dtc[num_rows != 1] - - # error if only fish with one observation. - if (nrow(dtc) == 0) stop("must have two observations to interpolate") - - # extract and determine start time - t_seq <- seq( - start_time, max(dtc$detection_timestamp_utc), - int_time_stamp - ) - - # bin data by time interval and add bin to dtc - dtc[, bin := t_seq[findInterval(detection_timestamp_utc, t_seq)]] - - # make all combinations of animals and detection bins - dtc <- dtc[data.table::CJ(bin = t_seq, animal_id = unique(animal_id)), - on = c("bin", "animal_id") - ] - - data.table::setkey(dtc, animal_id, bin, detection_timestamp_utc) - - # if only need to do linear interpolation: - if (is.null(trans) | lnl_thresh == 0) { - dtc[, bin_stamp := detection_timestamp_utc][ - is.na(detection_timestamp_utc), - bin_stamp := bin - ] - - dtc[, i_lat := approx(detection_timestamp_utc, deploy_lat, - xout = bin_stamp, ties = "ordered" - )$y, by = animal_id] - dtc[, i_lon := approx(detection_timestamp_utc, deploy_long, - xout = bin_stamp, ties = "ordered" - )$y, by = animal_id] - - dtc[is.na(deploy_long), record_type := "interpolated"] - dtc <- dtc[, c("animal_id", "bin_stamp", "i_lat", "i_lon", "record_type")] - det <- det[num_rows == 1, c( - "animal_id", "bin_stamp", "i_lat", "i_lon", - "record_type" - )] - out <- data.table::rbindlist(list(dtc, det), use.names = TRUE) - data.table::setkey(out, animal_id, bin_stamp) - out[, bin_stamp := t_seq[findInterval(bin_stamp, t_seq)]] - out <- na.omit(out, cols = "i_lat") - data.table::setnames(out, c( - "animal_id", "bin_timestamp", "latitude", - "longitude", "record_type" - )) - out <- unique(out) - out <- data.table::setorder(out, animal_id, bin_timestamp, -record_type) - - # If out_class == NULL, then return data as data.table - if (is.null(out_class)) { - out <- as.data.frame(out) - return(out) - } - - # if out_class == "tibble", then return tibble object - if (out_class == "tibble") { - out <- tibble::as_tibble(out) - return(out) - } - - # if out_class == NULL, then return data.frame object - return(out) - } - - # routine for combined nln and ln interpolation - - # identify start and end rows for observations before and after NA - ends <- dtc[!is.na(deploy_lat), .(start = .I[-nrow(.SD)], end = .I[-1]), - by = animal_id - ][end - start > 1] - - # identify observations that are both start and ends - dups <- c(ends$start, ends$end)[ends[, duplicated(c(start, end))]] - - # create and append duplicate rows for observations - # that are both start and end. - # This is so each observation can be in only one group - - # identifies rows and duplicate rows that need duplicated - dtc[, c("rep", "num") := list(1L, 1:.N)][dups, rep := 2L] - dtc <- dtc[rep(num, rep)] - dtc[, rep := NULL] - dtc[, num := NULL] - - # recalculate first and last rows- no duplicate rows this time... - new_ends <- dtc[!is.na(deploy_lat), .(start = .I[-nrow(.SD)], end = .I[-1]), - by = animal_id - ][end - start > 1] - - # create row index - dtc[, start_dtc := 1:.N] - - # extract rows that need interpolated - dtc <- dtc[new_ends, .( - animal_id = x.animal_id, - detection_timestamp_utc = x.detection_timestamp_utc, - deploy_lat = x.deploy_lat, deploy_long = x.deploy_long, - record_type = x.record_type, num_rows = x.num_rows, - bin = x.bin, i.start = start - ), - on = .(start_dtc >= start, start_dtc <= end) - ] - - # calculate great circle distance between coords - # dtc[, gcd := geosphere::distHaversine(as.matrix( - # .SD[1, c("deploy_long", "deploy_lat")]), - # as.matrix(.SD[.N, c("deploy_long", "deploy_lat")])), by = i.start] - - # remove geosphere dependency and use geodist. geosphere relies on sp and raster. - dtc[, gcd := geodist::geodist_vec(x1 = .SD[[1, "deploy_long"]], - y1 = .SD[[1, "deploy_lat"]], - x2 = .SD[[.N, "deploy_long"]], - y2 = .SD[[.N, "deploy_lat"]], - measure = "haversine"), - by = i.start] - - # calculate least cost (non-linear) distance between points - message("Calculating least-cost (non-linear) distances... (step 1 of 3)") - - grpn <- data.table::uniqueN(dtc$i.start) - - if (show_progress) pb <- txtProgressBar(min = 0, max = grpn, style = 3) - - dtc[, lcd := { - if (show_progress) setTxtProgressBar(pb, value = .GRP) - gdistance::costDistance(trans, - fromCoords = as.matrix( - .SD[1, c("deploy_long", "deploy_lat")] - ), - toCoords = as.matrix(.SD[.N, c("deploy_long", "deploy_lat")]) - ) - }, - by = i.start - ] - - # calculate ratio of gcd:lcd - dtc[, crit := gcd / lcd] - - # create keys for lookup - dtc[!is.na(detection_timestamp_utc), - t_lat := data.table::shift(deploy_lat, type = "lead"), - by = i.start - ] - dtc[!is.na(detection_timestamp_utc), - t_lon := data.table::shift(deploy_long, type = "lead"), - by = i.start - ] - dtc[!is.na(detection_timestamp_utc), - t_timestamp := data.table::shift(detection_timestamp_utc, type = "lead"), - by = i.start - ] - - # extract rows that need non-linear interpolation - # based on gcd:lcd distance - nln <- dtc[crit < lnl_thresh] - - land_chk <- dtc[is.infinite(lcd)][ - !is.na(deploy_lat), - c("deploy_lat", "deploy_long") - ] - - # stop execution and display offending receivers if any receivers are on land. - - capture <- function(x) paste(capture.output(print(x)), collapse = "\n") - - if (nrow(land_chk) > 0) { - stop("Some coordinates are on land or beyond extent. - Interpolation impossible! Check receiver locations or extents of transition - layer:\n", capture(as.data.table(land_chk)), call. = FALSE) - } - - # extract data for linear interpolation - # check to make sure that all points to be interpolated - # are within the tranition layer is needed before any interpolation. - - ln <- dtc[crit >= lnl_thresh | is.nan(crit)] - if (nrow(ln) == 0) { - ln <- data.table::data.table( - animal_id = character(), i_lat = numeric(), - i_lon = numeric(), - bin_stamp = as.POSIXct(character()), - record_type = character() - ) - } else { - message("Starting linear interpolation... (step 2 of 3)") - # linear interpolation - grpn <- uniqueN(ln$i.start) - if (show_progress) pb <- txtProgressBar(min = 0, max = grpn, style = 3) - ln[, bin_stamp := detection_timestamp_utc][ - is.na(detection_timestamp_utc), - bin_stamp := bin - ] - ln[, i_lat := { - if (show_progress) setTxtProgressBar(pb, .GRP) - tmp <- .SD[ - c(1, .N), - c("detection_timestamp_utc", "deploy_lat") - ] - approx(c(tmp$detection_timestamp_utc), - c(tmp$deploy_lat), - xout = c(bin_stamp), - ties = "ordered" - )$y - }, - by = i.start - ] - ln[, i_lon := { - tmp <- .SD[ - c(1, .N), - c("detection_timestamp_utc", "deploy_long") - ] - approx(c(tmp$detection_timestamp_utc), - c(tmp$deploy_long), - xout = c(bin_stamp), - ties = "ordered" - )$y - }, - by = i.start - ] - ln[is.na(deploy_long), record_type := "interpolated"] - } - - # extract records to lookup - nln_small <- nln[!is.na(detection_timestamp_utc)][!is.na(t_lat)] - - if (nrow(nln_small) == 0) { - nln <- data.table( - animal_id = character(), i_lat = numeric(), - i_lon = numeric(), - bin_stamp = as.POSIXct(character()), - record_type = character() - ) - } else { - # nln interpolation - # create lookup table - data.table::setkey(nln_small, deploy_lat, deploy_long, t_lat, t_lon) - lookup <- unique(nln_small[, .(deploy_lat, deploy_long, t_lat, t_lon), - allow.cartesian = TRUE - ]) - - message("\nStarting non-linear interpolation... (step 3 of 3)") - grpn <- nrow(lookup) - if (show_progress) pb <- txtProgressBar(min = 0, max = grpn, style = 3) - # calculate non-linear interpolation for all unique movements in lookup (old version with sp::coordinates) - ## lookup[, coord := { if(show_progress) setTxtProgressBar(pb, value = .GRP); - ## sp::coordinates( - ## gdistance::shortestPath(trans, as.matrix( - ## .SD[, c("deploy_long", "deploy_lat")]), as.matrix( - ## .SD[, c("t_lon", "t_lat")]), output = "SpatialLines") - ## ) - ## }, - ## by = 1:nrow(lookup)] - - # alternative using sf (needs tested!) - # st_as_sf(gdistance::shortestPath(trans, as.matrix(lookup[, c("deploy_long", "deploy_lat")]), as.matrix(lookup[, c("t_lon", "t_lat")]), output = "SpatialLines"), crs = 4326) - - lookup[, coord := { - if (show_progress) setTxtProgressBar(pb, value = .GRP) - sf::st_as_sf( - gdistance::shortestPath(trans, as.matrix( - .SD[, c("deploy_long", "deploy_lat")] - ), as.matrix( - .SD[, c("t_lon", "t_lat")] - ), output = "SpatialLines"), - crs = 4326 - ) - }, - by = 1:nrow(lookup) - ] - - message("\nFinalizing results.") - - lookup[, grp := 1:.N] - - # extract interpolated points from coordinate lists... - res <- lookup[, .( - nln_longitude = lookup$coord[[.I]][, 1], - nln_latitude = lookup$coord[[.I]][, 2] - ), by = grp] - - # set keys, join interpolation and original data - data.table::setkey(lookup, grp) - data.table::setkey(res, grp) - lookup <- lookup[res] - lookup[, coord := NULL] - - # added first/last rows, number sequence for groups - lookup[lookup[, .I[1], by = grp]$V1, nln_longitude := deploy_long] - lookup[lookup[, .I[.N], by = grp]$V1, nln_longitude := t_lon] - lookup[lookup[, .I[1], by = grp]$V1, nln_latitude := deploy_lat] - lookup[lookup[, .I[.N], by = grp]$V1, nln_latitude := t_lat] - lookup[, seq_count := 1:.N, by = grp] - - # lookup interpolated values for original dataset - data.table::setkey(lookup, deploy_lat, deploy_long, t_lat, t_lon) - nln_small <- lookup[nln_small, allow.cartesian = TRUE] - data.table::setkey(nln_small, i.start, seq_count) - - # add timeseries for interpolating nln movements - nln_small[ - nln_small[, .I[1], by = i.start]$V1, - i_time := detection_timestamp_utc - ] - nln_small[nln_small[, .I[.N], by = i.start]$V1, i_time := t_timestamp] - - # arch <- nln_small - # nln_small <- nln_small[i.start == 163] - - nln_small[, latitude_lead := data.table::shift(nln_latitude, type = "lag", fill = NA), by = i.start] - nln_small[, longitude_lead := data.table::shift(nln_longitude, type = "lag", fill = NA), by = i.start] - - # Switch from geosphere to geodist - # nln_small[, cumdist := geosphere::distGeo( - # .SD[, c("nln_longitude", "nln_latitude")], - # .SD[,c("longitude_lead", "latitude_lead")]), by = i.start] - - nln_small[, cumdist := geodist::geodist_vec(x1 = .SD[["nln_longitude"]], - y1 = .SD[["nln_latitude"]], - x2 = .SD[["longitude_lead"]], - y2 = .SD[["latitude_lead"]], - measure = 'geodesic'), - by = .I] - - nln_small[is.na(cumdist), cumdist := 0] - nln_small[, cumdist := cumsum(cumdist), by = i.start] - - nln_small[, latitude_lead := NULL][, longitude_lead := NULL] - - # calculate cumdist - ## nln_small[, cumdist := cumsum(c(0, sqrt(diff(nln_longitude) ^ 2 + - ## diff(nln_latitude) ^ 2))), - ## by = i.start] - - # interpolate missing timestamps for interpolated coordinates - nln_small[, i_time := as.POSIXct( - approx(cumdist, i_time, - xout = cumdist, - ties = "ordered" - )$y, - origin = "1970-01-01 00:00:00", - tz = attr(nln_small$i_time, "tzone") - ), - by = i.start - ] - - # create timestamp vector to interpolate on. - nln[, bin_stamp := detection_timestamp_utc] - nln[is.na(detection_timestamp_utc), bin_stamp := bin] - nln[, grp := i.start] - - # interpolate timestamps - data.table::setkey(nln_small, i.start) - data.table::setkey(nln, i.start) - nln[, i_lat := { - tmp <- nln_small[ - .(.SD[1, "i.start"]), - c("i_time", "nln_latitude") - ] - approx(tmp$i_time, tmp$nln_latitude, - xout = bin_stamp, - ties = "ordered" - )$y - }, - by = grp - ] - - nln[, i_lon := { - tmp <- nln_small[ - .(.SD[1, "i.start"]), - c("i_time", "nln_longitude") - ] - approx(tmp$i_time, tmp$nln_longitude, - xout = bin_stamp, - ties = "ordered" - )$y - }, - by = grp - ] - - nln[is.na(deploy_long), record_type := "interpolated"] - } - - # combine into a single data.table - out <- data.table::rbindlist(list( - ln[ - record_type == "interpolated", - c("animal_id", "bin_stamp", "i_lat", "i_lon", "record_type") - ], - nln[ - record_type == "interpolated", - c( - "animal_id", "bin_stamp", "i_lat", "i_lon", - "record_type" - ) - ], - det[, c( - "animal_id", "bin_stamp", "i_lat", "i_lon", - "record_type" - )] - ), use.names = TRUE) - - out[, !c("animal_id")] - data.table::setkey(out, animal_id, bin_stamp) - out[, bin_stamp := t_seq[findInterval(bin_stamp, t_seq)]] - data.table::setnames(out, c( - "animal_id", "bin_timestamp", "latitude", - "longitude", "record_type" - )) - out <- na.omit(out, cols = "latitude") - out <- unique(out) - data.table::setorder(out, animal_id, bin_timestamp, -record_type) - - # If out_class == NULL, then return data as data.table - if (is.null(out_class)) { - out <- as.data.frame(out) - return(out) - } - - # if out_class == "tibble", then return tibble object - if (out_class == "tibble") { - out <- tibble::as_tibble(out) - return(out) - } - - # if out_class == NULL, then return data.frame object - return(out) -} +#' Interpolate new positions within a spatiotemporal path data +#' +#' Interpolate new positions within a spatiotemporal path data set +#' (e.g., detections of tagged fish) at regularly-spaced time intervals +#' using linear or non-linear interpolation. +#' +#' @param det An object of class `glatos_detections` or data frame +#' containing spatiotemporal data with at least 4 columns containing +#' 'animal_id', 'detection_timestamp_utc', 'deploy_long', and +#' 'deploy_lat' columns. +#' +#' @param start_time specify the first time bin for interpolated data. +#' If not supplied, default is first timestamp in the input data +#' set. Must be a character string that can be coerced to +#' 'POSIXct' or an object of class 'POSIXct'. If character string +#' is supplied, timezone is automatically set to UTC. +#' +#' @param out_class Return results as a data.table or tibble. Default +#' returns results as data.frame. Accepts `data.table` or `tibble`. +#' +#' @param int_time_stamp The time step size (in seconds) of interpolated +#' positions. Default is 86400 (one day). +#' +#' @param trans An optional transition matrix with the "cost" of +#' moving across each cell within the map extent. Must be of class +#' `TransitionLayer`. A transition layer may be +#' created from a polygon shapefile using [make_transition]. +#' +#' @param lnl_thresh A numeric threshold for determining if linear or +#' non-linear interpolation shortest path will be used. +#' +#' @param show_progress Logical. Progress bar and status messages will be +#' shown if TRUE (default) and not shown if FALSE. +#' +#' @details Non-linear interpolation uses the `gdistance` package +#' to find the shortest pathway between two locations (i.e., +#' receivers) that avoid 'impossible' movements (e.g., over land for +#' fish). The shortest non-linear path between two locations is +#' calculated using a transition matrix layer that represents the +#' 'cost' of an animal moving between adjacent grid cells. The +#' transition matrix layer (see [gdistance]) is created from +#' a polygon shapefile using [make_transition] or from a +#' `RasterLayer` object using [transition][gdistance::transition]. In +#' `make_transition`, each cell in the output transition layer +#' is coded as water (1) or land (0) to represent possible (1) and +#' impossible (0) movement paths. +#' +#' @details Linear interpolation is used for all points when +#' `trans` is not supplied. When `trans` is supplied, +#' then interpolation method is determined for each pair of +#' sequential observed detections. For example, linear interpolation +#' will be used if the two geographical positions are exactly the +#' same and when the ratio (linear distance:non-linear distance) +#' between two positions is less than `lnl_thresh`. Non-linear +#' interpolation will be used when ratio is greater than +#' `lnl_thresh`. When the ratio of linear distance to +#' non-linear distance is greater than `lnl_thresh`, then the +#' distance of the non-linear path needed to avoid land is greater +#' than the linear path that crosses land. `lnl_thresh` can be +#' used to control whether non-linear or linear interpolation is +#' used for all points. For example, non-linear interpolation will +#' be used for all points when `lnl_thresh` > 1 and linear +#' interpolation will be used for all points when `lnl_thresh` +#' = 0. +#' +#' @details All linear interpolation is done by code{stats::approx} with +#' argument `ties = "ordered"` controlling how tied `x` values +#' are handled. See [approxfun()]. +#' +#' @return A dataframe with animal_id, bin_timestamp, +#' latitude, longitude, and record_type. +#' +#' +#' @author Todd Hayden, Tom Binder, Chris Holbrook +#' +#' @examples +#' +#' #-------------------------------------------------- +#' # EXAMPLE #1 - simple interpolate among lakes +#' +#' # get polygon of the Great Lakes +#' data(great_lakes_polygon) # glatos example data +#' plot(sf::st_geometry(great_lakes_polygon), xlim = c(-92, -76)) +#' +#' # make sample detections data frame +#' pos <- data.frame( +#' animal_id = 1, +#' deploy_long = c(-87, -82.5, -78), +#' deploy_lat = c(44, 44.5, 43.5), +#' detection_timestamp_utc = as.POSIXct(c( +#' "2000-01-01 00:00", +#' "2000-02-01 00:00", "2000-03-01 00:00" +#' ), tz = "UTC") +#' ) +#' +#' # add to plot +#' points(deploy_lat ~ deploy_long, data = pos, pch = 20, cex = 2, col = "red") +#' +#' # interpolate path using linear method +#' path1 <- interpolate_path(pos) +#' nrow(path1) # now 61 points +#' sum(path1$record_type == "interpolated") # 58 interpolated points +#' +#' # add linear path to plot +#' points(latitude ~ longitude, data = path1, pch = 20, cex = 0.8, col = "blue") +#' +#' # load a transition matrix of Great Lakes +#' # NOTE: This is a LOW RESOLUTION TransitionLayer suitable only for +#' # coarse/large scale interpolation only. Most realistic uses +#' # will need to create a TransitionLayer; see ?make_transition. +#' data(greatLakesTrLayer) # glatos example data; a TransitionLayer +#' +#' # interpolate path using non-linear method (requires 'trans') +#' path2 <- interpolate_path(pos, trans = greatLakesTrLayer) +#' +#' # add non-linear path to plot +#' points(latitude ~ longitude, +#' data = path2, pch = 20, cex = 1, +#' col = "green" +#' ) +#' +#' # can also force linear-interpolation with lnlThresh = 0 +#' path3 <- interpolate_path(pos, trans = greatLakesTrLayer, lnl_thresh = 0) +#' +#' # add new linear path to plot +#' points(latitude ~ longitude, +#' data = path3, pch = 20, cex = 1, +#' col = "magenta" +#' ) +#' +#' #-------------------------------------------------- +#' # EXAMPLE #2 - walleye in western Lake Erie +#' \dontrun{ +#' +#' +#' # get example walleye detection data +#' det_file <- system.file("extdata", "walleye_detections.csv", +#' package = "glatos" +#' ) +#' det <- read_glatos_detections(det_file) +#' +#' # take a look +#' head(det) +#' +#' # extract one fish and subset date +#' det <- det[det$animal_id == 22 & +#' det$detection_timestamp_utc > as.POSIXct("2012-04-08") & +#' det$detection_timestamp_utc < as.POSIXct("2013-04-15"), ] +#' +#' # get polygon of the Great Lakes +#' data(great_lakes_polygon) # glatos example data; an sf object +#' +#' # convert polygon to terra::spatVector +#' great_lakes_polygon <- terra::vect(great_lakes_polygon) +#' +#' # crop polygon to western Lake Erie +#' maumee <- terra::crop(great_lakes_polygon, +#' y = terra::ext(-83.7, -82.5, 41.3, 42.4) +#' ) +#' +#' +#' plot(maumee, col = "grey") +#' points(deploy_lat ~ deploy_long, +#' data = det, pch = 20, col = "red", +#' xlim = c(-83.7, -80) +#' ) +#' +#' +#' # make transition layer object +#' tran <- make_transition3(sf::st_as_sf(maumee), res = c(0.1, 0.1)) +#' +#' # plot to check output +#' plot(tran$rast, xlim = c(-83.7, -82.0), ylim = c(41.3, 42.7)) +#' plot(maumee, add = TRUE) +#' +#' # not high enough resolution- bump up resolution, will take some time +#' tran1 <- make_transition3(sf::st_as_sf(maumee), res = c(0.001, 0.001)) +#' +#' # plot to check resolution- much better +#' plot(tran1$rast, xlim = c(-83.7, -82.0), ylim = c(41.3, 42.7)) +#' plot(maumee, add = TRUE) +#' +#' +#' # add fish detections to make sure they are "on the map" +#' # plot unique values only for simplicity +#' foo <- unique(det[, c("deploy_lat", "deploy_long")]) +#' points(foo$deploy_long, foo$deploy_lat, pch = 20, col = "red") +#' +#' # call with "transition matrix" (non-linear interpolation), other options +#' # note that it is quite a bit slower than linear interpolation +#' pos2 <- interpolate_path(det, +#' trans = tran1$transition, +#' out_class = "data.table" +#' ) +#' +#' plot(maumee, col = "grey") +#' points(latitude ~ longitude, data = pos2, pch = 20, col = "red", cex = 0.5) +#' } +#' +#' @export + + + +interpolate_path <- function(det, trans = NULL, start_time = NULL, + int_time_stamp = 86400, lnl_thresh = 0.9, + out_class = NULL, show_progress = TRUE) { + # stop if out_class is not NULL, data.table, or tibble + if (!is.null(out_class)) { + if (!(out_class %in% c("data.table", "tibble"))) { + stop('out_class is not a "data.table" or "tibble"') + } + } + + # check to see that trans is a transition layer or transition stack + if (!is.null(trans) & + inherits(trans, c("TransitionLayer", "TransitionStack")) == FALSE) { + stop( + paste0( + "Supplied object for 'trans' argument is not class ", + "TransitionLayer or TransitionStack." + ), + call. = FALSE + ) + } + + # check start_time + if (is.null(start_time)) { + start_time <- min(det$detection_timestamp_utc) + } + if (is.na(start_time) & length(start_time) > 0) { + stop("start_time cannot be coerced to 'POSIXct' or 'POSIXt' class") + } + + if (is.character(start_time)) { + start_time <- as.POSIXct(start_time, tz = "UTC") + } + + # make sure start_time < largest timestamp in dataset + if (start_time > max(det$detection_timestamp_utc)) { + stop("start_time is larger than last detection. No data to interpolate!", call. = FALSE) + } + + # make copy of detections for function + dtc <- data.table::as.data.table(det) + + # subset only columns for function and rows >= start_time: + dtc <- dtc[detection_timestamp_utc >= start_time, c( + "animal_id", + "detection_timestamp_utc", + "deploy_lat", + "deploy_long" + )] + + dtc[, record_type := "detection"] + + # count number of rows- single observations are not interpolated + dtc[, num_rows := nrow(.SD), by = animal_id] + + # Sort detections by transmitter id and then by detection timestamp + data.table::setkey(dtc, animal_id, detection_timestamp_utc) + + # save original dataset to combine with interpolated data in the end + det <- data.table::copy(dtc) + data.table::setnames(det, c( + "animal_id", "bin_stamp", "i_lat", "i_lon", + "record_type", "num_rows" + )) + + # remove any fish with only one detection + dtc <- dtc[num_rows != 1] + + # error if only fish with one observation. + if (nrow(dtc) == 0) stop("must have two observations to interpolate") + + # extract and determine start time + t_seq <- seq( + start_time, max(dtc$detection_timestamp_utc), + int_time_stamp + ) + + # bin data by time interval and add bin to dtc + dtc[, bin := t_seq[findInterval(detection_timestamp_utc, t_seq)]] + + # make all combinations of animals and detection bins + dtc <- dtc[data.table::CJ(bin = t_seq, animal_id = unique(animal_id)), + on = c("bin", "animal_id") + ] + + data.table::setkey(dtc, animal_id, bin, detection_timestamp_utc) + + # if only need to do linear interpolation: + if (is.null(trans) | lnl_thresh == 0) { + dtc[, bin_stamp := detection_timestamp_utc][ + is.na(detection_timestamp_utc), + bin_stamp := bin + ] + + dtc[, i_lat := approx(detection_timestamp_utc, deploy_lat, + xout = bin_stamp, ties = "ordered" + )$y, by = animal_id] + dtc[, i_lon := approx(detection_timestamp_utc, deploy_long, + xout = bin_stamp, ties = "ordered" + )$y, by = animal_id] + + dtc[is.na(deploy_long), record_type := "interpolated"] + dtc <- dtc[, c("animal_id", "bin_stamp", "i_lat", "i_lon", "record_type")] + det <- det[num_rows == 1, c( + "animal_id", "bin_stamp", "i_lat", "i_lon", + "record_type" + )] + out <- data.table::rbindlist(list(dtc, det), use.names = TRUE) + data.table::setkey(out, animal_id, bin_stamp) + out[, bin_stamp := t_seq[findInterval(bin_stamp, t_seq)]] + out <- na.omit(out, cols = "i_lat") + data.table::setnames(out, c( + "animal_id", "bin_timestamp", "latitude", + "longitude", "record_type" + )) + out <- unique(out) + out <- data.table::setorder(out, animal_id, bin_timestamp, -record_type) + + # If out_class == NULL, then return data as data.table + if (is.null(out_class)) { + out <- as.data.frame(out) + return(out) + } + + # if out_class == "tibble", then return tibble object + if (out_class == "tibble") { + out <- tibble::as_tibble(out) + return(out) + } + + # if out_class == NULL, then return data.frame object + return(out) + } + + # routine for combined nln and ln interpolation + + # identify start and end rows for observations before and after NA + ends <- dtc[!is.na(deploy_lat), .(start = .I[-nrow(.SD)], end = .I[-1]), + by = animal_id + ][end - start > 1] + + # identify observations that are both start and ends + dups <- c(ends$start, ends$end)[ends[, duplicated(c(start, end))]] + + # create and append duplicate rows for observations + # that are both start and end. + # This is so each observation can be in only one group + + # identifies rows and duplicate rows that need duplicated + dtc[, c("rep", "num") := list(1L, 1:.N)][dups, rep := 2L] + dtc <- dtc[rep(num, rep)] + dtc[, rep := NULL] + dtc[, num := NULL] + + # recalculate first and last rows- no duplicate rows this time... + new_ends <- dtc[!is.na(deploy_lat), .(start = .I[-nrow(.SD)], end = .I[-1]), + by = animal_id + ][end - start > 1] + + # create row index + dtc[, start_dtc := 1:.N] + + # extract rows that need interpolated + dtc <- dtc[new_ends, .( + animal_id = x.animal_id, + detection_timestamp_utc = x.detection_timestamp_utc, + deploy_lat = x.deploy_lat, deploy_long = x.deploy_long, + record_type = x.record_type, num_rows = x.num_rows, + bin = x.bin, i.start = start + ), + on = .(start_dtc >= start, start_dtc <= end) + ] + + # calculate great circle distance between coords + # dtc[, gcd := geosphere::distHaversine(as.matrix( + # .SD[1, c("deploy_long", "deploy_lat")]), + # as.matrix(.SD[.N, c("deploy_long", "deploy_lat")])), by = i.start] + + # remove geosphere dependency and use geodist. geosphere relies on sp and raster. + dtc[, gcd := geodist::geodist_vec( + x1 = .SD[[1, "deploy_long"]], + y1 = .SD[[1, "deploy_lat"]], + x2 = .SD[[.N, "deploy_long"]], + y2 = .SD[[.N, "deploy_lat"]], + measure = "haversine" + ), + by = i.start + ] + + # calculate least cost (non-linear) distance between points + message("Calculating least-cost (non-linear) distances... (step 1 of 3)") + + grpn <- data.table::uniqueN(dtc$i.start) + + if (show_progress) pb <- txtProgressBar(min = 0, max = grpn, style = 3) + + dtc[, lcd := { + if (show_progress) setTxtProgressBar(pb, value = .GRP) + gdistance::costDistance(trans, + fromCoords = as.matrix( + .SD[1, c("deploy_long", "deploy_lat")] + ), + toCoords = as.matrix(.SD[.N, c("deploy_long", "deploy_lat")]) + ) + }, + by = i.start + ] + + # calculate ratio of gcd:lcd + dtc[, crit := gcd / lcd] + + # create keys for lookup + dtc[!is.na(detection_timestamp_utc), + t_lat := data.table::shift(deploy_lat, type = "lead"), + by = i.start + ] + dtc[!is.na(detection_timestamp_utc), + t_lon := data.table::shift(deploy_long, type = "lead"), + by = i.start + ] + dtc[!is.na(detection_timestamp_utc), + t_timestamp := data.table::shift(detection_timestamp_utc, type = "lead"), + by = i.start + ] + + # extract rows that need non-linear interpolation + # based on gcd:lcd distance + nln <- dtc[crit < lnl_thresh] + + land_chk <- dtc[is.infinite(lcd)][ + !is.na(deploy_lat), + c("deploy_lat", "deploy_long") + ] + + # stop execution and display offending receivers if any receivers are on land. + + capture <- function(x) paste(capture.output(print(x)), collapse = "\n") + + if (nrow(land_chk) > 0) { + stop("Some coordinates are on land or beyond extent. + Interpolation impossible! Check receiver locations or extents of transition + layer:\n", capture(as.data.table(land_chk)), call. = FALSE) + } + + # extract data for linear interpolation + # check to make sure that all points to be interpolated + # are within the tranition layer is needed before any interpolation. + + ln <- dtc[crit >= lnl_thresh | is.nan(crit)] + if (nrow(ln) == 0) { + ln <- data.table::data.table( + animal_id = character(), i_lat = numeric(), + i_lon = numeric(), + bin_stamp = as.POSIXct(character()), + record_type = character() + ) + } else { + message("Starting linear interpolation... (step 2 of 3)") + # linear interpolation + grpn <- uniqueN(ln$i.start) + if (show_progress) pb <- txtProgressBar(min = 0, max = grpn, style = 3) + ln[, bin_stamp := detection_timestamp_utc][ + is.na(detection_timestamp_utc), + bin_stamp := bin + ] + ln[, i_lat := { + if (show_progress) setTxtProgressBar(pb, .GRP) + tmp <- .SD[ + c(1, .N), + c("detection_timestamp_utc", "deploy_lat") + ] + approx(c(tmp$detection_timestamp_utc), + c(tmp$deploy_lat), + xout = c(bin_stamp), + ties = "ordered" + )$y + }, + by = i.start + ] + ln[, i_lon := { + tmp <- .SD[ + c(1, .N), + c("detection_timestamp_utc", "deploy_long") + ] + approx(c(tmp$detection_timestamp_utc), + c(tmp$deploy_long), + xout = c(bin_stamp), + ties = "ordered" + )$y + }, + by = i.start + ] + ln[is.na(deploy_long), record_type := "interpolated"] + } + + # extract records to lookup + nln_small <- nln[!is.na(detection_timestamp_utc)][!is.na(t_lat)] + + if (nrow(nln_small) == 0) { + nln <- data.table( + animal_id = character(), i_lat = numeric(), + i_lon = numeric(), + bin_stamp = as.POSIXct(character()), + record_type = character() + ) + } else { + # nln interpolation + # create lookup table + data.table::setkey(nln_small, deploy_lat, deploy_long, t_lat, t_lon) + lookup <- unique(nln_small[, .(deploy_lat, deploy_long, t_lat, t_lon), + allow.cartesian = TRUE + ]) + + message("\nStarting non-linear interpolation... (step 3 of 3)") + grpn <- nrow(lookup) + if (show_progress) pb <- txtProgressBar(min = 0, max = grpn, style = 3) + # calculate non-linear interpolation for all unique movements in lookup (old version with sp::coordinates) + ## lookup[, coord := { if(show_progress) setTxtProgressBar(pb, value = .GRP); + ## sp::coordinates( + ## gdistance::shortestPath(trans, as.matrix( + ## .SD[, c("deploy_long", "deploy_lat")]), as.matrix( + ## .SD[, c("t_lon", "t_lat")]), output = "SpatialLines") + ## ) + ## }, + ## by = 1:nrow(lookup)] + + # alternative using sf (needs tested!) + # st_as_sf(gdistance::shortestPath(trans, as.matrix(lookup[, c("deploy_long", "deploy_lat")]), as.matrix(lookup[, c("t_lon", "t_lat")]), output = "SpatialLines"), crs = 4326) + + lookup[, coord := { + if (show_progress) setTxtProgressBar(pb, value = .GRP) + sf::st_as_sf( + gdistance::shortestPath(trans, as.matrix( + .SD[, c("deploy_long", "deploy_lat")] + ), as.matrix( + .SD[, c("t_lon", "t_lat")] + ), output = "SpatialLines"), + crs = 4326 + ) + }, + by = 1:nrow(lookup) + ] + + message("\nFinalizing results.") + + lookup[, grp := 1:.N] + + # extract interpolated points from coordinate lists... + res <- lookup[, .( + nln_longitude = lookup$coord[[.I]][, 1], + nln_latitude = lookup$coord[[.I]][, 2] + ), by = grp] + + # set keys, join interpolation and original data + data.table::setkey(lookup, grp) + data.table::setkey(res, grp) + lookup <- lookup[res] + lookup[, coord := NULL] + + # added first/last rows, number sequence for groups + lookup[lookup[, .I[1], by = grp]$V1, nln_longitude := deploy_long] + lookup[lookup[, .I[.N], by = grp]$V1, nln_longitude := t_lon] + lookup[lookup[, .I[1], by = grp]$V1, nln_latitude := deploy_lat] + lookup[lookup[, .I[.N], by = grp]$V1, nln_latitude := t_lat] + lookup[, seq_count := 1:.N, by = grp] + + # lookup interpolated values for original dataset + data.table::setkey(lookup, deploy_lat, deploy_long, t_lat, t_lon) + nln_small <- lookup[nln_small, allow.cartesian = TRUE] + data.table::setkey(nln_small, i.start, seq_count) + + # add timeseries for interpolating nln movements + nln_small[ + nln_small[, .I[1], by = i.start]$V1, + i_time := detection_timestamp_utc + ] + nln_small[nln_small[, .I[.N], by = i.start]$V1, i_time := t_timestamp] + + # arch <- nln_small + # nln_small <- nln_small[i.start == 163] + + nln_small[, latitude_lead := data.table::shift(nln_latitude, type = "lag", fill = NA), by = i.start] + nln_small[, longitude_lead := data.table::shift(nln_longitude, type = "lag", fill = NA), by = i.start] + + # Switch from geosphere to geodist + # nln_small[, cumdist := geosphere::distGeo( + # .SD[, c("nln_longitude", "nln_latitude")], + # .SD[,c("longitude_lead", "latitude_lead")]), by = i.start] + + nln_small[, cumdist := geodist::geodist_vec( + x1 = .SD[["nln_longitude"]], + y1 = .SD[["nln_latitude"]], + x2 = .SD[["longitude_lead"]], + y2 = .SD[["latitude_lead"]], + measure = "geodesic" + ), + by = .I + ] + + nln_small[is.na(cumdist), cumdist := 0] + nln_small[, cumdist := cumsum(cumdist), by = i.start] + + nln_small[, latitude_lead := NULL][, longitude_lead := NULL] + + # calculate cumdist + ## nln_small[, cumdist := cumsum(c(0, sqrt(diff(nln_longitude) ^ 2 + + ## diff(nln_latitude) ^ 2))), + ## by = i.start] + + # interpolate missing timestamps for interpolated coordinates + nln_small[, i_time := as.POSIXct( + approx(cumdist, i_time, + xout = cumdist, + ties = "ordered" + )$y, + origin = "1970-01-01 00:00:00", + tz = attr(nln_small$i_time, "tzone") + ), + by = i.start + ] + + # create timestamp vector to interpolate on. + nln[, bin_stamp := detection_timestamp_utc] + nln[is.na(detection_timestamp_utc), bin_stamp := bin] + nln[, grp := i.start] + + # interpolate timestamps + data.table::setkey(nln_small, i.start) + data.table::setkey(nln, i.start) + nln[, i_lat := { + tmp <- nln_small[ + .(.SD[1, "i.start"]), + c("i_time", "nln_latitude") + ] + approx(tmp$i_time, tmp$nln_latitude, + xout = bin_stamp, + ties = "ordered" + )$y + }, + by = grp + ] + + nln[, i_lon := { + tmp <- nln_small[ + .(.SD[1, "i.start"]), + c("i_time", "nln_longitude") + ] + approx(tmp$i_time, tmp$nln_longitude, + xout = bin_stamp, + ties = "ordered" + )$y + }, + by = grp + ] + + nln[is.na(deploy_long), record_type := "interpolated"] + } + + # combine into a single data.table + out <- data.table::rbindlist(list( + ln[ + record_type == "interpolated", + c("animal_id", "bin_stamp", "i_lat", "i_lon", "record_type") + ], + nln[ + record_type == "interpolated", + c( + "animal_id", "bin_stamp", "i_lat", "i_lon", + "record_type" + ) + ], + det[, c( + "animal_id", "bin_stamp", "i_lat", "i_lon", + "record_type" + )] + ), use.names = TRUE) + + out[, !c("animal_id")] + data.table::setkey(out, animal_id, bin_stamp) + out[, bin_stamp := t_seq[findInterval(bin_stamp, t_seq)]] + data.table::setnames(out, c( + "animal_id", "bin_timestamp", "latitude", + "longitude", "record_type" + )) + out <- na.omit(out, cols = "latitude") + out <- unique(out) + data.table::setorder(out, animal_id, bin_timestamp, -record_type) + + # If out_class == NULL, then return data as data.table + if (is.null(out_class)) { + out <- as.data.frame(out) + return(out) + } + + # if out_class == "tibble", then return tibble object + if (out_class == "tibble") { + out <- tibble::as_tibble(out) + return(out) + } + + # if out_class == NULL, then return data.frame object + return(out) +}