R/helper_functions.R

In khufkens/swiftr: Process bird (swift) GPS position tracker data and meta-data.

#' Converts altitudes in m to hPa
#'
#' uses sealevel mean values as a reference
#'
#' Reference / source needs verification:
#' https://www.engineeringtoolbox.com/air-altitude-pressure-d_462.html
#'
#' @param m vector of altitudes in meter
#'
#' @return altitude expressed as hPa
#' @export

meter_to_hPa <- function(m){
  if(missing(m)){
    stop("please provide an altitude")
  }
  (101325 * (1 - (2.25577 * 10^-5 * m)) ^ 5.25588) / 100
}

#' Calculate wind direction from CDS data
#'
#' @param nc CDS netcdf data file
#'
#' @return absolute wind speed and direction
#' @export
#'
#' @examples
#'

wind_direction <- function(nc){

  # read in data
  u_ms <- raster(nc, varname = "u")
  v_ms <- raster(nc, varname = "v")

  # calculate absolute windspeed
  wind_speed <- sqrt(u_ms^2 + v_ms^2)

  # in radians
  wind_direction <- atan2(u_ms, v_ms)  * 180/pi

  # regularize to 0 - 360
  wind_direction <- (wind_direction + 360) %% 360

  return(stack(wind_speed, wind_direction))
}


# Calculates the geodesic distance and bearing between two points
# specified by degrees latitude/longitude using the Haversine formula (hf)
# https://en.wikipedia.org/wiki/Haversine_formula

great_circle_distance <- function(lat1, lon1, lat2, lon2){

  # convert to radian
  lat1 <- lat1 * pi/180
  lat2 <- lat2 * pi/180
  lon1 <- lon1 * pi/180
  lon2 <- lon2 * pi/180

  # Earth mean radius [km]
  R <- 6371

  # calculate delta values
  d_lon <- (lon2 - lon1)
  d_lat <- (lat2 - lat1)

  # haversine
  a <- sin(d_lat/2)^2 + cos(lat1) * cos(lat2) * sin(d_lon/2)^2
  c <- 2 * asin(min(1,sqrt(a)))
  d <- R * c

  # convert to m
  d <- d * 1000

  x <- sin(d_lon) * cos(lat2)
  y <- cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(d_lon)

  # calculate heading
  heading <- atan2(y, x) * 180/pi
  heading <- (heading + 360) %% 360

  return(data.frame(distance = d,
                    heading = heading))
}