mass updates 4 hi-res acoustic telemetry, etc

.Rbuildignore

@@ -26,3 +26,4 @@
 ^CRAN-SUBMISSION$
 ^doc$
 ^Meta$
+^data-raw$

DESCRIPTION

@@ -29,7 +29,7 @@ License: MIT + file LICENSE
 LazyData: true
 NeedsCompilation: yes
 Encoding: UTF-8
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Imports: 
     tibble (>= 2.1.3),
     ggplot2 (>= 3.4.0),
@@ -41,16 +41,15 @@ Imports:
     CircStats,
     mvtnorm,
     tmvtnorm,
-    raster,
-    broom,
     rnaturalearth,
     traipse,
     pathroutr (>= 0.2.1),
     rnaturalearthhires (>= 0.2.0),
     rnaturalearthdata,
     rosm,
     ggspatial,
-    mapproj
+    mapproj,
+    terra
 LinkingTo: TMB, RcppEigen
 Suggests: 
     testthat,

NAMESPACE

@@ -31,7 +31,6 @@ importFrom(TMB,MakeADFun)
 importFrom(TMB,newtonOption)
 importFrom(TMB,oneStepPredict)
 importFrom(TMB,sdreport)
-importFrom(broom,tidy)
 importFrom(dplyr,"%>%")
 importFrom(dplyr,arrange)
 importFrom(dplyr,bind_cols)
@@ -57,8 +56,6 @@ importFrom(dplyr,ungroup)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,aes_string)
 importFrom(ggplot2,coord_fixed)
-importFrom(ggplot2,coord_map)
-importFrom(ggplot2,coord_quickmap)
 importFrom(ggplot2,coord_sf)
 importFrom(ggplot2,element_blank)
 importFrom(ggplot2,element_line)
@@ -72,7 +69,6 @@ importFrom(ggplot2,geom_hline)
 importFrom(ggplot2,geom_line)
 importFrom(ggplot2,geom_path)
 importFrom(ggplot2,geom_point)
-importFrom(ggplot2,geom_polygon)
 importFrom(ggplot2,geom_qq)
 importFrom(ggplot2,geom_qq_line)
 importFrom(ggplot2,geom_ribbon)
@@ -97,7 +93,6 @@ importFrom(ggplot2,scale_size)
 importFrom(ggplot2,theme)
 importFrom(ggplot2,theme_bw)
 importFrom(ggplot2,theme_minimal)
-importFrom(ggplot2,theme_void)
 importFrom(ggplot2,unit)
 importFrom(ggplot2,vars)
 importFrom(ggplot2,xlab)
@@ -111,9 +106,6 @@ importFrom(grDevices,hcl.colors)
 importFrom(grDevices,hcl.pals)
 importFrom(mvtnorm,rmvnorm)
 importFrom(patchwork,wrap_plots)
-importFrom(raster,extent)
-importFrom(raster,extract)
-importFrom(raster,nlayers)
 importFrom(rnaturalearth,ne_countries)
 importFrom(sf,"st_geometry<-")
 importFrom(sf,st_as_sf)
@@ -143,6 +135,7 @@ importFrom(stats,BIC)
 importFrom(stats,acf)
 importFrom(stats,approx)
 importFrom(stats,cov)
+importFrom(stats,filter)
 importFrom(stats,median)
 importFrom(stats,na.omit)
 importFrom(stats,nlminb)
@@ -158,6 +151,8 @@ importFrom(stats,rgamma)
 importFrom(stats,rlnorm)
 importFrom(stats,runif)
 importFrom(stats,sd)
+importFrom(terra,ext)
+importFrom(terra,extract)
 importFrom(tibble,as_tibble)
 importFrom(tibble,tibble)
 importFrom(tidyr,gather)

R/aniMotum-package.R

@@ -29,7 +29,7 @@
 ##' @importFrom tibble as_tibble tibble
 ##' @importFrom sf st_as_sf st_set_crs st_transform st_is_longlat st_crs
 ##' @importFrom sf st_coordinates st_geometry<- st_bbox st_cast
-##' @importFrom raster extent extract nlayers
+##' @importFrom terra ext extract
 ##' @importFrom traipse track_distance track_angle track_distance_to
 ##' @importFrom TMB MakeADFun sdreport newtonOption oneStepPredict
 ##' @importFrom stats approx cov sd predict nlminb optim na.omit median qlogis qnorm pnorm runif
@@ -39,7 +39,6 @@
 ##' @importFrom ggplot2 element_text scale_colour_manual scale_colour_gradientn
 ##' @importFrom grDevices extendrange grey hcl.colors
 ##' @importFrom rnaturalearth ne_countries
-##' @importFrom broom tidy
 NULL
 
 ##' @name ellie

R/fit_ssm.R

@@ -4,19 +4,20 @@
 ##' @description fits: i) a simple random walk (`rw`) ii) a correlated random walk
 ##' (`crw` - a random walk on velocity), or iii) a time-varying move persistence 
 ##' model (`mp`), all in continuous-time, to filter Argos LS, and/or KF/KS 
-##' location data, processed light-level geolocation data (GLS), and/or GPS data. 
-##' Location data of different types can combined in a single data frame 
-##' (see details). Predicts locations at user-specified time intervals 
-##' (regular or irregular).
+##' location data, GPS data, and/or generic locations with associated standard 
+##' errors (e.g., processed light-level geolocation data, or high-resolution 
+##' acoustic telemetry data). Location data of different types can combined in a 
+##' single data frame (see details). Predicts locations at user-specified time 
+##' intervals (regular or irregular).
 ##'
 ##' @param x a `data.frame`, `tibble` or `sf-tibble` of observations, depending 
 ##' on the tracking data type. See more in the Details section, below, and the 
 ##' Overview vignette \code{vignette("Overview", package = "aniMotum")}.
-##' @param vmax max travel rate (m/s) passed to [trip::sda] to identify
+##' @param vmax max travel rate (m/s) to identify
 ##'  outlier locations
 ##' @param ang angles (deg) of outlier location "spikes" 
 ##' @param distlim lengths (m) of outlier location "spikes" 
-##' @param spdf (logical) turn [trip::sda] pre-filtering on (default; TRUE) or off
+##' @param spdf (logical) turn pre-filtering on (default; TRUE) or off
 ##' @param min.dt minimum allowable time difference between observations;
 ##' `dt <= min.dt` will be ignored by the SSM
 ##' @param pf just pre-filter the data, do not fit the SSM (default is FALSE)
@@ -51,10 +52,11 @@
 ##' string in YYYY-MM-DD HH:MM:SS format. If a text string is supplied then the
 ##' time zone is assumed to be `UTC`. lc (location class) can include the 
 ##' following values: 3, 2, 1, 0, A, B, Z, G, or GL. The latter two are for GPS 
-##' and GLS locations, respectively. Class Z values are assumed to have the same 
-##' error variances as class B. By default, class `G` (GPS) locations are assumed 
-##' to have error variances 10x smaller than Argos class 3 variances, but unlike 
-##' Argos error variances the GPS variances are the same for longitude and latitude. 
+##' locations and 'Generic Locations', respectively. Class Z values are assumed 
+##' to have the same error variances as class B. By default, class `G` (GPS) 
+##' locations are assumed to have error variances 10x smaller than Argos class 3 
+##' variances, but unlike Argos error variances the GPS variances are the same for 
+##' longitude and latitude. 
 ##' 
 ##' The [aniMotum::format_data] function can be used as a data pre-processing 
 ##' step or called automatically within `fit_ssm` to restructure data that is 
@@ -69,10 +71,14 @@
 ##' including the above 5 plus **`smaj`, `smin`, `eor`** that contain Argos error 
 ##' ellipse variables (in m for `smaj`, `smin` and deg for `eor`). 
 ##' 
-##' Light-level geolocation (GLS) locations can be modelled provided each 
-##' longitude and latitude has a corresponding standard error. These data should 
-##' have 7 columns, including the above 5 plus `lonerr`, `laterr` (in degrees). 
-##' In this case, all lc values should be set to `GL`. 
+##' Generic locations can be modelled provided each longitude and latitude 
+##' (or X and Y) coordinate has a corresponding standard error. These data should 
+##' have 7 columns, including the above 5 plus two extra columns, typically
+##' named `x.sd`, `y.sd` that provide the standard errors for the longitude, 
+##' latitude (or X, Y) coordinates. Longitude and latitude standard errors should
+##' be in degrees, whereas X and Y standard errors should be in m. In either case, 
+##' all `lc` values should be set to `GL` (Generic Location), the helper function
+##' [format_data] will add the `lc` variable to the input data automatically. 
 ##' 
 ##' Multiple location data types can be combined in a single data frame 
 ##' (see the Overview vignette for examples). 

R/format_data.R

@@ -24,15 +24,16 @@
 ##' "smin" (ellipse semi-minor axis), and "eor" (ellipse orientation). Ignored if
 ##' these variables are missing from the input data.
 ##' @param sderr the names (as quoted character strings) of provided standard 
-##' errors in longitude and latitude: defaults are "lonerr", "laterr". Typically,
-##' these are only provided for processed light-level geolocation data. Ignored if
-##' these variables are missing from the input data.
+##' errors for `lon,lat` or `x,y`: default names are `x.sd`, `y.sd`. Typically,
+##' these are only provided for generic location data such as processed light-level
+##' geolocations, or high-resolution acoustic detections. The argument is ignored
+##' if these variables are missing from the input data.
 ##' @param tz the timezone the applies to the data/time variable if they are not 
 ##' in `tz = 'UTC'`. A list of valid timezone names can be viewed via `OlsonNames()`
 ##' 
 ##' @return a data.frame or sf-tibble of input data in expected aniMotum format. 
 ##' Additional columns required by `fit_ssm()`, if missing, will be added to the 
-##' formatted tibble: `smaj`, `smin`, `eor`, `lonerr`, and `laterr`.
+##' formatted tibble: `smaj`, `smin`, `eor`, `x.sd`, and `y.sd`.
 ##' 
 ##' @importFrom sf st_crs
 ##' @importFrom dplyr tibble select everything
@@ -57,7 +58,7 @@ format_data <- function(x,
                         lc = "lc",
                         coord = c("lon","lat"),
                         epar = c("smaj","smin","eor"),
-                        sderr = c("lonerr","laterr"),
+                        sderr = c("x.sd","y.sd"),
                         tz = "UTC") {
 
   ## check that all variable names are character strings
@@ -100,7 +101,7 @@ format_data <- function(x,
 
   ## add lc if missing from input data
   if (!lc %in% names(x)) {
-    ## Case when data are GLS/geolocations
+    ## Case when data are Generic Locations
     if (all(!epar %in% names(x)) & all(sderr %in% names(x))) {
       if (inherits(x, "data.frame", which = TRUE) == 1) {
         x <- data.frame(x, lc = rep("GL", nrow(x)))
@@ -169,18 +170,18 @@ format_data <- function(x,
     }
   }
   if(all(!epar %in% names(x), sderr %in% names(x))) {
-    ## GLS data
+    ## Generic Location data
     ## add expected error ellipse variables
     x$smaj <- x$smin <- x$eor <- as.double(NA)
     xx <- x[, c(id, date, lc, coord, epar, sderr, xt.vars)]
     if(!inherits(x, "sf")) {
-      names(xx)[c(1:5, 9:10)] <- c("id","date","lc",coord,"lonerr","laterr")
-      names(xx)[4:5] <- c("lon","lat")
+      names(xx)[c(1:5, 9:10)] <- c("id","date","lc",coord,"x.sd","y.sd")
+      #names(xx)[4:5] <- c("lon","lat")
     } else if (inherits(x, "sf")) {
-      names(xx)[c(1:4, 8:9)] <- c("id","date","lc",coord,"lonerr","laterr")
+      names(xx)[c(1:4, 8:9)] <- c("id","date","lc",coord,"x.sd","y.sd")
     }
   }
-  
+
   ## in cases where user supplies id as a factor, drop any unused factor levels 
   ##    and coerce to character
   if(is.factor(xx$id)) xx$id <- droplevels(xx$id)

R/mpfilter.R

@@ -28,7 +28,7 @@
 ##' @param lpsi is deprecated, use ssm_control(lower = list(lpsi = -Inf)) instead, see \code{ssm_control} for details
 ##'
 ##' @importFrom TMB MakeADFun sdreport newtonOption FreeADFun
-##' @importFrom stats approx cov sd predict nlminb optim na.omit
+##' @importFrom stats approx cov sd predict nlminb optim na.omit filter
 ##' @importFrom utils flush.console
 ##' @importFrom tibble as_tibble
 ##' @importFrom sf st_crs st_coordinates st_geometry<- st_as_sf st_set_crs
@@ -145,7 +145,7 @@ mpfilter <-
              xout = d.all$date,
              rule = 2)$y
     x.init <-
-      as.numeric(stats::filter(x.init1, rep(1, 5) / 5))
+      as.numeric(filter(x.init1, rep(1, 5) / 5))
     x.na <- which(is.na(x.init))
     x.init[x.na] <- x.init1[x.na]
 
@@ -185,8 +185,8 @@ mpfilter <-
 
     ## start to work out which obs_mod to use for each observation
     d <- d %>% mutate(obs.type = factor(obs.type, 
-                                        levels = c("LS","KF","GPS"), 
-                                        labels = c("LS","KF","GPS"))
+                                        levels = c("LS","KF","GL","GPS"), 
+                                        labels = c("LS","KF","GL","GPS"))
                       )
     p.obst <- table(d$obs.type) / nrow(d)
     obst <- which(table(d$obs.type) > 0)
@@ -244,7 +244,8 @@ mpfilter <-
       m = d.all$smin,
       M = d.all$smaj,
       c = d.all$eor,
-      K = cbind(d.all$emf.x, d.all$emf.y)
+      K = cbind(d.all$emf.x, d.all$emf.y),
+      GLerr = cbind(d.all$x.sd, d.all$y.sd)
     ) 
 
     ## TMB - create objective function

R/pf_add_emf.R

@@ -22,12 +22,18 @@ pf_add_emf <- function(x, emf) {
   x$lc <- with(x, as.character(lc))
   x <- merge(x, tmp, by = "lc", all.x = TRUE, sort = FALSE)
 
+  if(all("lonerr" %in% names(x), "laterr" %in% names(x))) {
+    x <- x[order(x$date), c("id","date","lc","smaj","smin","eor",
+                            "lonerr","laterr","keep","obs.type",
+                            "emf.x","emf.y","geometry")]  
+  } else {
+    x <- x[order(x$date), c("id","date","lc","smaj","smin","eor",
+                            "x.sd","y.sd","keep","obs.type",
+                            "emf.x","emf.y","geometry")]
+  }
 
-  x <- x[order(x$date), c("id","date","lc","smaj","smin","eor",
-                                "lonerr","laterr","keep","obs.type",
-                                "emf.x","emf.y","geometry")]
-  x$emf.x <- with(x, ifelse(obs.type %in% c("KF","GLS"), NA, emf.x))
-  x$emf.y <- with(x, ifelse(obs.type %in% c("KF","GLS"), NA, emf.y))
+  x$emf.x <- with(x, ifelse(obs.type %in% c("KF","GL"), NA, emf.x))
+  x$emf.y <- with(x, ifelse(obs.type %in% c("KF","GL"), NA, emf.y))
 
   if (sum(is.na(x$lc)) > 0)
     stop(

R/pf_dup_dates.R

@@ -9,10 +9,9 @@
 pf_dup_dates <- function(x, min.dt) {
 
   ##  flag any duplicate date records,
-
-  x$keep <- with(x, difftime(date, c(as.POSIXct(NA), date[-nrow(x)]), 
-                             units = "secs") > min.dt)
-  x$keep <- with(x, ifelse(is.na(keep), TRUE, keep))
+  x$keep <- difftime(x$date, c(as.POSIXct(NA), x$date[-nrow(x)]), 
+                     units = "secs") > min.dt
+  x$keep <- ifelse(is.na(x$keep), TRUE, x$keep)
 
   return(x)
 }

R/pf_obs_type.R

@@ -19,22 +19,30 @@ pf_obs_type <- function(x) {
                                "GPS", obs.type))
   x$obs.type <- with(x, ifelse(lc == "GL" & 
                                  (is.na(smaj) | is.na(smin) | is.na(eor)) & 
-                                 (!is.na(lonerr) & !is.na(laterr)), 
-                               "GLS", obs.type))
-  
+                                 (!is.na(x.sd) & !is.na(y.sd)), 
+                               "GL", obs.type))
+
   ##  if any records with smaj/smin = 0 then set to NA and obs.type to "LS"
   ##  convert error ellipse smaj & smin from m to km and eor from deg to rad
   x$smaj <- with(x, ifelse(smaj == 0 | smin == 0, NA, smaj)) / 1000
   x$smin <- with(x, ifelse(smin == 0 | is.na(smaj), NA, smin)) / 1000
   x$eor <- with(x, ifelse(is.na(smaj) & is.na(smin), NA, eor)) / 180 * pi
 
   x$obs.type <- with(x, ifelse(is.na(smaj) & is.na(smin) & is.na(eor) & 
-                                 (obs.type != "GLS" & obs.type != "GPS"), 
+                                 (obs.type != "GL" & obs.type != "GPS"), 
                                "LS", obs.type))
 
-  ## convert GLS errors from degrees lon/lat to km
-  x$lonerr <- with(x, lonerr * 6378.137 / 180 * pi)  
-  x$laterr <- with(x, laterr * 6378.137 / 180 * pi)
+  if(all("lon" %in% names(x), "lat" %in% names(x), "lonerr" %in% names(x), "laterr" %in% names(x))) {
+    ## if GL SD's are loneer/laterr then convert from deg to km
+    x$x.sd <- with(x, x.sd * 6378.137 / 180 * pi)  
+    x$y.sd <- with(x, y.sd * 6378.137 / 180 * pi)
+
+  } else if(all((all("lon" %in% names(x), "lat" %in% names(x)) |
+             all("x" %in% names(x), "y" %in% names(x))), "x.sd" %in% names(x), "y.sd" %in% names(x))){
+    ## if GL SD's are x.sd/y.sd then convert from m to km
+    x$x.sd <- with(x, x.sd / 1000)  
+    x$y.sd <- with(x, y.sd / 1000)
+  }
 
   return(x)
 }