earthtools: Tools for working on planet Earth

#' compute location of coordinates after great circle projection
#' 
#' This function provides a convienant wrapper to 
#' \code{\link[geosphere]{destPoint}}
#' 
#' 
#' @param x An object used to determine which implementation to use
#' @param ... It's an S3 thing. You wouldn't understand.
#' @param .data An object that inherits from \code{\link[base]{data.frame}}. In 
#'   general this will be on of \code{data.frame}, 
#'   \code{\link[data.table]{data.table}}, or \code{\link[dplyr]{tbl_df}}
#' @param latitude Either a numeric vector of latitudes [degrees] or the column
#'   of \code{.data} which contains latitudes. This maybe quoted or unquoted;
#'   see examples.
#' @param longitude Either a numeric vector of longitudes [degrees]
#'   corresponding with latitude or the column of \code{.data} which contains
#'   longitudes. This maybe quoted or unquoted; see examples.
#' @param bearing Either a numeric vector of bearings [degrees] or the column of
#'   \code{.data} which contains bearings/headings. This maybe quoted or
#'   unquoted see examples.
#' @param distance Either a numeric vector of projection distance [nautical
#'   miles] or the column of \code{.data} which contains projection distances.
#'   This maybe quoted or unquoted see examples.
#' @param output_type string in \code{c("matrix", "data.table", "data.frame", 
#'   "list")}
#' @param method Either \code{"GC"} [default] or \code{"rhumb"}. Used to declare
#'   either a great circle calculation or rhumb line calculation
#' @return If \code{.data} is supplied, an object of the same type and with the
#'   same columns as \code{.data} plus two more, \code{end_latitude} and
#'   \code{end_longitude}.  Otherwise, an object of type determined by
#'   output_type which will generally have two columns, latitude and longitude.
#'   If the input coordinates have length 1, then a named numeric vector is
#'   returned.
#'   
#' @examples
#' # basic use
#' compute_projection(39.86167, -104.6732, 90, 15)
#' compute_projection(39.86167, -104.6732, 86:90, 1:15)
#' 
#' # use inside a data.table
#' library(data.table)
#' apts <- data.table(airport=c("ATL", "DEN", "ORD", "SEA"),
#'                    latitude=c(33.63670, 39.86167, 41.97933, 47.44989),
#'                    longitude=c(-84.42786, -104.67317, -87.90739, -122.31178))
#' apts[, c("platitude", "plongitude"):=compute_projection(latitude, longitude, 90, 15)]
#' 
#' # use with magrittr
#' library(magrittr)
#' apts %>% compute_projection(latitude, longitude, 90, 15)
#' 
#' # columns as strings
#' lat_col <- names(apts)[2]
#' apts %>% compute_projection(lat_col, "longitude", 90, 15)
#' 
#' # predict next position
#' tracks <- data.frame(id = c("a","b","c"),
#'                      lat = 0,
#'                      lon = 0,
#'                      heading = 30,
#'                      ground_speed = seq(300,360, 30))
#' time_step <- 1/60 #one minute
#' 
#' tracks %>% compute_projection(lat, lon, heading, tracks$ground_speed*time_step)
#' 
#' @export
#' @importFrom geosphere destPoint
#' @rdname compute_projection
compute_projection <- function(x, ..., method="GC"){
  UseMethod("compute_projection")
}

#' @export
#' @rdname compute_projection
#' @importFrom dplyr rename bind_cols %>%
compute_projection.data.frame <- function(.data, latitude, longitude, bearing, distance, method="GC"){
  lat_ <- determine_val(.data, substitute(latitude))
  lon_ <- determine_val(.data, substitute(longitude))
  bearing_ <- determine_val(.data, substitute(bearing))
  dist_ <- determine_val(.data, substitute(distance))
  
  compute_projection.numeric(lat_, lon_, bearing_, dist_,output_type = class(.data), method=method) %>%
    rename(end_latitude=latitude, end_longitude=longitude) %>%
    bind_cols(.data,.)%>%
    format_return(class(.data)) %>%
    return
}

# attempt to return vectors from .data or the value provided by supplied argument
determine_val <- function(.data, val){
  cls <- class(val)
  if(cls == "call"){
    val <- eval(val)
    cls <- class(val)
  }  
  if(cls == "name"){
    val <- ifelse(as.character(val) %in% names(.data), as.character(val), eval(val))
  }
  if(cls %in% c("integer","numeric") || is.numeric(val)){
    return(val)
  }
  return(.data[[which(names(.data)==val)]])
}

format_return <- function(.data, return_type){
  if("data.table" %in% return_type){
    data.table::setDT(.data)
    return(.data)
  }
  return(as(object = .data, Class = return_type[1]))
  
}

#' @export
#' @importFrom data.table data.table
#' @importFrom dplyr tbl tbl_df
#' @importFrom geosphere destPoint destPointRhumb
#' @rdname compute_projection
compute_projection.numeric <- function(latitude, longitude, bearing, distance, output_type="data.table", method="GC") {
  valid_out_types <- c("matrix", "data.table", "data.frame", "list", "tbl", "tbl_df")
  if (!any(output_type %in% valid_out_types)) {
    stop(paste0("invalid output type specified, must be one of:\n\t", 
                paste0(valid_out_types, collapse = ", ")))
  }
  
  if(method=="GC")
   points <- destPoint(p=matrix(c(longitude, latitude), ncol = 2, byrow = FALSE), 
                      b=bearing, d=distance, r=.earth.radius.NMi)
  else if (method == "rhumb")
   points <- destPointRhumb(p=matrix(c(longitude, latitude), ncol = 2, byrow = FALSE), 
                                       b=bearing, d=distance, r=.earth.radius.NMi)
  else
    stop("unknown method: ", method, ". Use either GC or rhumb")
  
  if ("matrix" %in% output_type) {
    colnames(points) <- c("longitude", "latitude")
    return(points[,2:1])
  } else if ("data.table" %in% output_type) {
    return(data.table(latitude=points[, "lat"], longitude=points[, "lon"]))
  } else if ("tbl_df" %in% output_type) {
    return(tbl_df(data.frame(latitude=points[, "lat"], longitude=points[, "lon"])))
  } else if ("tbl" %in% output_type) {
    return(tbl(data.frame(latitude=points[, "lat"], longitude=points[, "lon"])))
  } else if ("data.frame" %in% output_type) {
    return(data.frame(latitude=points[, "lat"], longitude=points[, "lon"]))
  } else if ("list" %in% output_type) {
    return(list(list(latitude=points[, "lat"]), list(longitude=points[, "lon"])))
  } 
}

#' compute location of coordinates after rhumb line projection
#' 
#' @description 
#' **DEPRECATED** Please use \code{compute_projection(method="rhumb")}
#' 
#' This function provides a convienant wrapper to 
#' \code{\link[geosphere]{destPointRhumb}}
#' 
#' @param latitude numeric vector of latitudes [degrees]
#' @param longitude numeric vector of longitudes [degrees] corresponding with latitude
#' @param bearing initial bearing [degrees]
#' @param distance projection distance [nautical miles]
#' @param output_type string in \code{c("matrix", "data.table", "data.frame", "list")}
#' @importFrom geosphere destPointRhumb
#' @return a \code{data.table} with columns latitude and longitude. If the input coordinates have length 1, then a named numeric vector is returned if 
#' @export
#' 
compute_projection_rhumb <- function(latitude, longitude, bearing, distance, output_type="data.table") {
  warning("compute_projection_rhumb is deprecated and will go away when earthtools 2.x.y is released.  Use compute_projection(method=\"rhumb\") instead")
 compute_projection.numeric(latitude, longitude, bearing, distance, output_type, method = "rhumb")
}

#' produce a ring of pseudo-nodes
#'
#' It can be useful to add a range ring of nodes to a graph. This function will select such
#' nodes subject to radius and disance constraints. As a general rule the max distance should
#' be twice the max fix distance tolerance so that aircraft "pass by" no more or less than one
#' fix on the ring (assuming they don't double back or fly in at a near tangent heading)
#'
#' @param center_latitude Numeric. The latitude of the range ring center point.
#' @param center_longitude Numeric. The longitude of the range ring center point.
#' @param radius Numeric. The radius of the arc about the center point in nautical miles.
#' @param dist_between_points Numeric. The maximum distance between points on the
#'   range ring. Actual distance is between this upper bound and this minus the tolerance,
#'   all in nautical miles.
#' @param tolerance Numeric. How close to the target distance between points do you need to get?
#' @param from_angle Numeric. The angle from the center point to the arc defining the left-most
#'   point of the arc.
#' @param to_angle Numeric. The angle from the center point to the arc defining the right-most
#'   point of the arc (moving clockwise from the \code{from_angle}).
#' 
#' @export
project_arc <- function(center_latitude, center_longitude, radius, dist_between_points, tolerance=0.01,
                        from_angle = 0, to_angle = 359.999) {
  
  theta0 <- asin(dist_between_points/radius) * 180/pi
  theta <- search_for_theta(interval=c(theta0-.5*theta0, theta0+.5*theta0),
                            expected_separation=dist_between_points,
                            tolerance=tolerance,
                            center_latitude=center_latitude,
                            center_longitude=center_longitude,
                            radius=radius)
  
  if (from_angle > to_angle)
    to_angle <- to_angle + 360
  
  return(compute_projection(center_latitude, center_longitude, 
                            bearing=seq(from_angle, to_angle, by=theta) %% 360,
                            distance=radius))
}

#' @noRd
#' 
#' Basically perform a binomial search over the interval to find an angle increment
#' that will space points along the arc (w.r.t. great circle distance) by the expected
#' separation
search_for_theta <- function(interval, expected_separation, tolerance,
                             center_latitude, center_longitude, radius) {
  theta <- mean(interval)
  p <- compute_projection(center_latitude, center_longitude, bearing=c(0, theta), distance=radius)
  d <- compute_dist_haversine(p[1, latitude], p[1, longitude],p[2, latitude], p[2, longitude])
  
  delta <- expected_separation-d
  if (delta > 0 && delta < tolerance)
    return(theta)
  
  if (delta < 0)
    return(search_for_theta(interval=c(min(interval), theta),
                            expected_separation=expected_separation,
                            tolerance=tolerance,
                            center_latitude=center_latitude,
                            center_longitude=center_longitude,
                            radius=radius))
  else # if (result < expected)
    return(search_for_theta(interval=c(theta, max(interval)),
                            expected_separation=expected_separation,
                            tolerance=tolerance,
                            center_latitude=center_latitude,
                            center_longitude=center_longitude,
                            radius=radius))
}