prefilter now calls component fn's;

add tz arg to format_data;
update Overview vignette;
sese2_n data dates changed to character;
increment version & date;

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: foieGras
 Title: Fit Continuous-Time State-Space and Latent Variable Models for Quality Control of Argos Satellite (and Other) Telemetry Data and for Estimating Movement Behaviour
-Version: 1.0-8.9128
-Date: 2022-11-08
+Version: 1.0-9
+Date: 2022-11-11
 Authors@R: 
     c(
       person(given = "Ian",

R/fit_ssm.R

@@ -67,7 +67,7 @@
 ##' In this case, all lc values should be set to `GL`. 
 ##' 
 ##' Multiple location data types can be combined in a single data frame 
-##' (see the vignette for examples). 
+##' (see the Overview vignette for examples). 
 ##' 
 ##' When data are provided as an `sf-tibble`, the user-specified projection is 
 ##' respected, although projected units are always transformed to km to improve 
@@ -81,7 +81,7 @@
 ##' * `predicted` an sf tbl of predicted location states
 ##' * `fitted` an sf tbl of fitted locations
 ##' * `par` model parameter summary
-##' * `data` an augmented sf tbl of the input data
+##' * `data` an augmented sf tbl of the formatted input data
 ##' * `inits` a list of initial values
 ##' * `pm` the process model fit, either "rw" or "crw"
 ##' * `ts` time time.step in h used

R/format_data.R

@@ -27,6 +27,8 @@
 ##' 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.
+##' @param tz the timezone the applies to the data/time variable if they are not 
+##' in `tz = 'UTC'`.
 ##' 
 ##' @return a data.frame or sf-tibble of input data in expected foieGras format. 
 ##' Additional columns required by `fit_ssm()`, if missing, will be added to the 
@@ -38,12 +40,13 @@
 ##' @examples
 ##' ## as a data pre-processing step
 ##' data(sese2_n)
-##' d <- format_data(sese2_n, date = "time", coord = c("longitude","latitude"))
+##' d <- format_data(sese2_n, date = "time", coord = c("longitude","latitude"), 
+##' tz = "America/Halifax")
 ##' fit <- fit_ssm(d, model = "crw", time.step = 24)
 ##' 
 ##' ## called automatically within fit_ssm()
-##' fit <- fit_ssm(sese2_n, date = "time", coord = c("longitude", "latitude"),
-##' model = "crw", time.step = 24)
+##' fit <- fit_ssm(sese2_n, date = "time", coord = c("longitude", "latitude"), 
+##' tz = "America/Halifax", model = "crw", time.step = 24)
 ##' @export
 ##' @md
 
@@ -53,7 +56,8 @@ format_data <- function(x,
                         lc = "lc",
                         coord = c("lon","lat"),
                         epar = c("smaj","smin","eor"),
-                        sderr = c("lonerr","laterr")) {
+                        sderr = c("lonerr","laterr"),
+                        tz = "UTC") {
 
   ## check that all variable names are character strings
   if(id %in% names(x)) 
@@ -171,10 +175,9 @@ format_data <- function(x,
   if(is.factor(xx$id)) xx$id <- droplevels(xx$id)
   xx$id <- as.character(xx$id)
 
-  ## convert dates to POSIXt if not already, assume tz = 'UTC'
+  ## convert dates to POSIXt if not already
   if(!inherits(xx$date, "POSIXt")) {
-    message("converting date variable to POSIX with timezone = 'UTC', if this timezone is incorrect then convert dates to POSIX format prior to using this function")
-    xx$date <- as.POSIXct(xx$date, tz = "UTC")
+    xx$date <- as.POSIXct(xx$date, tz = tz)
   }
   ## order records by date
   xx <- xx[order(xx$date), ]

R/pf_add_emf.R

@@ -0,0 +1,36 @@
+##' @title Add error multiplication factors for Argos LS data
+##' 
+##' @description Adds location error multiplication factors based on Argos location
+##' class (for type LS) & finalises prefiltered data for use by `sfilter()`
+##'
+##' @param x data from `pf_sf_project`
+##' @param emf optionally supplied data.frame of error multiplication factors 
+##' for Argos location quality classes. Relevant to Argos LS and GPS data only
+##' @keywords internal
+##' @md
+
+pf_add_emf <- function(x, emf) {
+
+  ## add LS error info to corresponding records
+  ## set emf's = NA if obs.type %in% c("KF","GL") - not essential but for clarity
+  if(is.null(emf)) {
+    tmp <- emf()
+  } else if(is.data.frame(emf)) {
+    tmp <- emf
+  } 
+
+  x$lc <- with(x, as.character(lc))
+  x <- merge(x, tmp, by = "lc", all.x = TRUE, sort = FALSE)
+  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))
+
+  if (sum(is.na(x$lc)) > 0)
+    stop(
+      "\n NA's found in location class values"
+    )
+
+  return(x)
+}

R/pf_dup_dates.R

@@ -0,0 +1,18 @@
+##' @title Find observations with duplicate dates
+##'
+##' @param x input data from `format_data()`
+##' @param min.dt minimum allowable time difference in s between observations; 
+##' `dt < min.dt` will be ignored by the SSM
+##' @keywords internal
+##' @md
+
+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))
+
+  return(x)
+}

R/pf_obs_type.R

@@ -0,0 +1,40 @@
+##' @title Determine observation type: LS, KF, GPS, or GLS
+##'
+##' @param x data from `pf_dup_dates()`
+##' @keywords internal
+##' @md
+
+pf_obs_type <- function(x) {
+
+  ## determine observation type: LS, KF, GPS or GLS
+  x$obs.type <- NA
+  x$obs.type <- with(x, 
+                     ifelse(!is.na(smaj) & !is.na(smin) & !is.na(eor), 
+                            "KF", obs.type))
+  x$obs.type <- with(x, ifelse(lc %in% c(3,2,1,0,"A","B","Z") & 
+                                 (is.na(smaj) | is.na(smin) | is.na(eor)), 
+                               "LS", obs.type))
+  x$obs.type <- with(x, ifelse(lc == "G" & 
+                                 (is.na(smaj) | is.na(smin) |is.na(eor)), 
+                               "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))
+
+  ##  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(any(is.na(smaj), is.na(smin), is.na(eor)) & 
+                                 all(obs.type != "GLS", 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)
+
+  return(x)
+}

R/pf_sda_filter.R

@@ -0,0 +1,96 @@
+##' @title Use `trip::sda` to identify extreme locations
+##' 
+##' @details `trip::sda` is a vectorized implementation of the Argos filter 
+##' presented in Freitas et al. (2008) Mar Mamm Sci 24:315-325. This function
+##' checks for errors returned by `trip::sda` and falls back to using the simpler 
+##' `trip::speedfilter` if an error is returned.
+##'
+##' @param x data from `pf_obs_type()`
+##' @param spdf turn speed filter on/off (logical; default is TRUE)
+##' @param vmax max travel rate (m/s)
+##' @param ang angles of outlier location "spikes" (default is `c(15,25)` deg);
+##' `ang = NA` turns off `trip::sda` filter in favour of 
+##' `trip::speedfilter`
+##' @param distlim lengths of outlier location "spikes" in km (default is 
+##' `c(2.5, 5)` m); `distlim = NA` turns off `trip::sda` filter 
+##' in favour of `trip::speedfilter`. Either `ang = NA` or 
+##' `distlim = NA` are sufficient.
+##' @importFrom trip sda speedfilter trip
+##' @importFrom sf st_coordinates st_is_longlat st_crs st_transform 
+##' @keywords internal
+##' @md
+
+pf_sda_filter <- function(x, spdf, vmax, ang, distlim) {
+
+## Use trip::sda to identify outlier locations
+if (spdf) {
+  if(inherits(x, "sf") && st_is_longlat(x)) {
+
+    xy <- as.data.frame(st_coordinates(x))
+    names(xy) <- c("lon","lat")
+    x <- cbind(x, xy)
+
+  } else if(inherits(x, "sf") && !st_is_longlat(x)) {
+
+    xy <- st_transform(x, crs = st_crs("+proj=longlat +datum=WGS84 +no_defs"))
+    xy <- as.data.frame(st_coordinates(xy))
+    names(xy) <- c("lon","lat")
+    x <- cbind(x, xy)
+  } 
+  x.tr <- subset(x, keep)[, c("lon","lat","date","id","lc","smaj","smin",
+                              "eor","lonerr","laterr","keep","obs.type")]
+  names(x.tr)[1:2] <- c("x","y")
+  x.tr <- suppressWarnings(trip(as.data.frame(x.tr), TORnames = c("date", "id"), 
+                                correct_all = FALSE))
+
+  if(any(is.na(ang))) ang <- c(0,0)
+  if(any(is.na(distlim))) distlim <- c(0,0)
+
+  tmp <-
+    suppressWarnings(try(
+      sda(
+        x.tr,
+        smax = vmax * 3.6,    # convert m/s to km/h
+        ang = ang,
+        distlim = distlim / 1000     # convert m to km
+      ),
+      silent = TRUE)
+    )
+  ## screen potential sdafilter errors
+  if (inherits(tmp, "try-error")) {
+
+    warning(
+      paste(
+        "\ntrip::sda produced an error on id",
+        x$id[1],
+        "using trip::speedfilter instead"
+      ),
+      immediate. = TRUE
+    )
+
+    tmp <-
+      suppressWarnings(try(
+        speedfilter(x.tr,
+                    max.speed = vmax * 3.6    # convert m/s to km/h
+        ),
+        silent = TRUE)
+      )
+
+    if (inherits(tmp, "try-error")) {
+
+      warning(
+        paste(
+          "\ntrip::speedfilter also produced an error on id",
+          x$id[1],
+          "can not apply speed filter prior to SSM filtering"
+        ),
+        immediate. = TRUE
+      )
+    }
+  }
+  x[x$keep, "keep"] <- tmp
+
+} 
+
+  return(x)
+}

R/pf_sf_project.R

@@ -0,0 +1,56 @@
+##' @title handle spatial projection
+##' 
+##' @description project from longlat to merc or respect user-supplied 
+##' projection (if not longlat) & ensure that longitudes straddling -180,180 or 
+##' 0,360 are shifted appropriately.
+##'
+##' @param x data from `pf_sda_filter`
+##' @importFrom sf st_as_sf st_crs st_transform st_is_longlat
+##' @keywords internal
+##' @md
+
+pf_sf_project <- function(x) {
+
+  if(!inherits(x, "sf")) {
+    ##  if lon spans -180,180 then shift to
+    ##    0,360; else if lon spans 360,0 then shift to
+    ##    -180,180 ... have to do this on keep subset only
+
+    xx <- subset(x, keep)
+
+    ## projection not provided by user so project to Mercator
+    sf_locs <- st_as_sf(x, coords = c("lon", "lat"), 
+                        crs = st_crs("+proj=longlat +datum=WGS84 +no_defs"))
+
+    if (any(diff(wrap_lon(xx$lon, 0)) > 300)) {
+      prj <- "+proj=merc +lon_0=0 +datum=WGS84 +units=km +no_defs"
+    } else if (any(diff(wrap_lon(xx$lon,-180)) < -300) ||
+               any(diff(wrap_lon(xx$lon,-180)) > 300)) {
+      prj <- "+proj=merc +lon_0=180 +datum=WGS84 +units=km +no_defs"
+    } else {
+      prj <- "+proj=merc +lon_0=0 +datum=WGS84 +units=km +no_defs"
+    }
+
+    sf_locs <-  st_transform(sf_locs, st_crs(prj))
+
+  } else {
+    ## if input data projection is longlat then set prj merc, otherwise respect 
+    ##     user-supplied projection
+    if(st_is_longlat(x)) {
+      prj <- "+proj=merc +lon_0=0 +datum=WGS84 +units=km +no_defs"
+    } else {
+      prj <- st_crs(x)$input  
+    }
+
+    # if data CRS units are m then change to km, otherwise optimiser may choke
+    if (grepl("units=m", prj, fixed = TRUE)) {
+      message("Converting projection units from m to km for efficient optimization")
+      prj <- sub("units=m", "units=km", prj, fixed = TRUE)
+    }
+    ll <- which(names(x) %in% c("lon","lat"))
+    sf_locs <- x[, -ll]
+    sf_locs <- st_transform(sf_locs, prj)
+  }
+
+  return(sf_locs)
+}