commit 06b1a6d31693fa73957be57d469d16d7fb9d414f
parent ed5b9602c10cfe5dcf4f024d2571a40c55d483d9
Author: eamoncaddigan <eamon.caddigan@gmail.com>
Date: Fri, 16 Sep 2016 15:24:27 -0400
Just copied the code over from flightconflicts, not integrated yet
Diffstat:
3 files changed, 155 insertions(+), 0 deletions(-)
diff --git a/R/createTrajectory.R b/R/createTrajectory.R
@@ -0,0 +1,90 @@
+#' Create a flighttrajectory object from the flight info.
+#'
+#' @param longitude Required; numeric vector giving aircraft longitude in
+#' degrees.
+#' @param latitude Required; numeric vector giving aircraft latitude in degrees.
+#' @param altitude Optional; numeric vector giving aircraft altitude (AGL) in
+#' feet. If missing, it will be set to 0.
+#' @param timestamp Optional; numeric vector giving the time of each observation
+#' in seconds. If missing, the observation period is assumed to be 1 s.
+#' @param bearing Optional; numeric vector giving the current bearing in
+#' degrees. If missing, it is estimated using pairs of successive lon/lat
+#' observations.
+#' @param groundspeed Optional; numeric vector giving the current ground speed
+#' of the aircraft in knots. If missing, it is estimated using pairs of
+#' successive lon/lat observations.
+#' @return A flighttrajectory object encapsulating these parameters (with
+#' default values substituded as necessary).
+#'
+#' @details \code{longitude} and \code{latitude} must be the same length.
+#' \code{timestamp}, \code{bearing}, and \code{groundspeed}, if present, must
+#' also match this length. \code{altitude} must also have a length equal to
+#' these parameters or be scalar.
+#'
+#' @export
+createTrajectory <- function(longitude, latitude, altitude = 0, timestamp = NULL,
+ bearing = NULL, groundspeed = NULL) {
+ if (!is.numeric(longitude)) stop("\"longitude\" must be a numeric vector")
+ nCoord <- length(longitude)
+
+ # Helper function to throw an error if the length of a vector is incorrect.
+ checkLength <- function(x) {
+ if (!is.numeric(x)) {
+ stop("\"", deparse(substitute(x)), "\" must be a numeric vector")
+ } else if (length(x) != nCoord) {
+ stop("Vector \"", deparse(substitute(x)), "\" has length = ", length(x),
+ ", expected length = ", nCoord)
+ }
+ return(TRUE)
+ }
+
+ checkLength(latitude)
+ coords <- cbind(longitude, latitude)
+
+ if (length(altitude) == 1) {
+ altitude <- rep(altitude, nCoord)
+ } else {
+ checkLength(altitude)
+ }
+
+ if (is.null(timestamp)) {
+ timestamp <- seq(1, nCoord)
+ } else {
+ checkLength(timestamp)
+ }
+
+ # Use flightpathr to calculate bearing between successive points if not
+ # specified.
+ if (is.null(bearing)) {
+ bearing <- flightpathr::coordsToBearing(cbind(coords, altitude))
+ bearing[nCoord] <- bearing[nCoord-1]
+ } else {
+ checkLength(bearing)
+ }
+
+ # Use geosphere to find the distance between points and use the timestamps to
+ # calculate groundspeed if not specified.
+ if (is.null(groundspeed)) {
+ distNM <- geosphere::distCosine(coords[1:(nCoord-1), ],
+ coords[2:nCoord, ],
+ r = 3444)
+ groundspeed <- distNM / diff(timestamp) * 3600
+ groundspeed <- c(groundspeed, groundspeed[nCoord-1])
+ } else{
+ checkLength(groundspeed)
+ }
+
+ flighttrajectory <- list(longitude = longitude,
+ latitude = latitude,
+ altitude = altitude,
+ timestamp = timestamp,
+ bearing = bearing,
+ groundspeed = groundspeed)
+ class(flighttrajectory) <- "flighttrajectory"
+
+ return(flighttrajectory)
+}
+
+#' Check if an object is a flighttrajectory
+#' @export
+is.flighttrajectory <- function(x) inherits(x, "flighttrajectory")
diff --git a/R/interpolateTrajectory.R b/R/interpolateTrajectory.R
@@ -0,0 +1,31 @@
+#' Interpolate a trajectory (in time)
+#'
+#' @param trajectory A \code{flighttrajectory} object.
+#' @param timestamp The new timestamp along which the data should be
+#' interpolated.
+#' @return A new \code{flighttrajectory} with the given \code{timestamp}
+#'
+#' @details This just performs linear interpolation for all of the values in the
+#' trajectory. A better approach would make use of the bearing and velocity
+#' information to smoothly interpolate the coordinates.
+#'
+#' @export
+interpolateTrajectory <- function(trajectory, timestamp) {
+ # We'll interpolate all of the features (except the timestamp, which is
+ # specified).
+ trajectoryFeatures <- names(trajectory)
+ trajectoryFeatures <- trajectoryFeatures[trajectoryFeatures != "timestamp"]
+
+ # Create a new trajectory with
+ newTrajectory <- list()
+ for (trajectoryFeature in trajectoryFeatures) {
+ newTrajectory[[trajectoryFeature]] <- approx(x = trajectory$timestamp,
+ y = trajectory[[trajectoryFeature]],
+ xout = timestamp,
+ method = "linear",
+ rule = 2)$y
+ }
+ newTrajectory[["timestamp"]] <- timestamp
+
+ return(do.call(createTrajectory, newTrajectory))
+}
diff --git a/tests/testthat/test_interpolateTrajectory.R b/tests/testthat/test_interpolateTrajectory.R
@@ -0,0 +1,34 @@
+library(flightconflicts)
+context("interpolateTrajectory")
+
+library(geosphere)
+
+kacy <- c(-74.5771667, 39.4575833)
+k17n <- c(-75.0330031, 39.7054758)
+
+# Two identical trajectories with different sampling rates
+coords1 <- gcIntermediate(kacy, k17n, n = 61)
+trajectory1 <-createTrajectory(coords1[, 1], coords1[, 2], altitude = 3500,
+ timestamp = seq(0, 800, length.out = 61))
+coords2 <- gcIntermediate(kacy, k17n, n = 229)
+trajectory2 <-createTrajectory(coords2[, 1], coords2[, 2], altitude = 3500,
+ timestamp = seq(0, 800, length.out = 229))
+
+test_that("Interpolated trajectory is close to correct", {
+ trajectoryInterpolated <- interpolateTrajectory(trajectory1,
+ trajectory2$timestamp)
+
+ for (tf in names(trajectoryInterpolated)) {
+ expect_equal(trajectoryInterpolated[[tf]], trajectory2[[tf]], tolerance = 1)
+ }
+})
+
+test_that("Correctly handling times outside original range", {
+ trajectoryInterpolated <- interpolateTrajectory(trajectory1,
+ c(-100, 0, 800, 900))
+
+ for (tf in names(trajectoryInterpolated)[names(trajectoryInterpolated) != "timestamp"]) {
+ expect_equal(trajectoryInterpolated[[tf]][1], trajectoryInterpolated[[tf]][2])
+ expect_equal(trajectoryInterpolated[[tf]][3], trajectoryInterpolated[[tf]][4])
+ }
+})
