library(readr)
library(dplyr)
library(sf)
library(ggplot2)
library(stringr)
date_of_interest <- "2017-08-01"
# get airport relevant details
source(here::here("R", "egll-data.R"))
arp <- get("egll_arp")
# use st_buffer to enlarge the "manually" defined holding polygons
# see https://stackoverflow.com/a/46706245/963575
library(units)
# buffer for holding polygon: 2 nautical miles
buf <- set_units(2, nmile) %>%
set_units(m)
holdings <- egll_holdings_sf %>%
st_transform(29902) %>%
st_buffer(buf) %>%
st_transform(4326)
# read flight info
flt_rt_file <- str_glue("data/egll_flt_rt_{d}.csv", d = date_of_interest)
cols <- cols(
flight_id = col_character(),
callsign = col_character(),
aircraft_reg = col_character(),
aircraft_type = col_character(),
aircraft_address = col_character(),
adep = col_character(),
ades = col_character(),
ssr_codes = col_character(),
period_start = col_datetime(format = ""),
period_finish = col_datetime(format = ""),
source_ids = col_character()
)
flt_rt <- read_csv(flt_rt_file)
flt_rt_no_dups <- flt_rt %>%
group_by(callsign) %>%
mutate(n = n()) %>%
ungroup() %>%
filter(n == 1 | (n >= 2 & !is.na(aircraft_type)))
# read positions within 40 NM
cols <- cols(
flight_id = col_character(),
timestamp = col_datetime(format = ""),
latitude = col_double(),
longitude = col_double(),
altitude = col_double(),
speed_gnd = col_double(),
track_gnd = col_double(),
vert_speed = col_double(),
on_ground = col_logical(),
distance = col_double(),
distance_arp = col_double()
)
pos_rt_40NM_file <- str_glue("data/egll_pos_rt_40NM_{d}.csv", d = date_of_interest)
pos_rt_40NM <- read_csv(pos_rt_40NM_file, col_types = cols) %>%
select(flight_id,
timestamp,
longitude, latitude, altitude,
distance, distance_arp) %>%
rename(distance_flown = distance) %>%
group_by(flight_id) %>%
arrange(timestamp) %>%
# add
# - sequence number
# - total of points
# convert
# - distance_arp from meters to NM
mutate(sequence_number = dplyr::row_number(),
n = n(),
distance_arp = distance_arp / 1852) %>%
ungroup()
#' extract reansition points at 40 NM, holdings and landing
#'
#' @param df A reference trajectory data frame
#' @param holdings_sf An sf polygon defining the various holding areas
#' @param min_fl Minimum flight level (with margin) for holdings.
#' For example at EGLL lowest flying level is at 7000 ft,
#' so in this function we set it at 6700...
#'
#' @return A data frame with points transitionsing from the various zones
extract_transition_points <- function(df, holdings_sf, min_fl = 6700) {
p <- df %>% group_by(flight_id)
# indexes within 40 NM
p_40 <- p %>%
filter(row_number() == 1) %>%
mutate(type = "P40") %>%
ungroup()
# indexes closest to ARP: an approximation of the landing spot
p_landing <- p %>%
filter(distance_arp == min(distance_arp)) %>%
filter(timestamp == min(timestamp)) %>%
mutate(type = "PLAND") %>%
ungroup()
# indexes for holding portions
st_as_xyz <- function(x) {
# TODO: stop if geometry is not XY POINT
data.frame(st_set_geometry(x, NULL), st_coordinates(x)) %>%
as_tibble() %>%
rename(longitude = X, latitude = Y)
}
# NOTE: for EGLL holding lower flight level is at 7000 ft
p_holding <- p %>%
st_as_sf(
coords = c("longitude", "latitude"),
dim = "XY",
crs = 4326) %>%
# NOTE: stay close minumum level in holding
filter(altitude > min_fl) %>%
st_intersection(holdings_sf) %>%
group_by(flight_id, id) %>%
filter(row_number() == 1 | row_number() == n()) %>%
mutate(type = "PHOLD") %>%
ungroup() %>%
st_as_xyz()
bind_rows(p_40, p_holding, p_landing) %>%
rename(holding_id = id)
}
transition_points <- extract_transition_points(pos_rt_40NM, holdings, min_fl = 6700)
# add bearing and quadrant
transition_points <- transition_points %>%
add_bearing(arp) %>%
add_quadrant()
# join positions with relevant flight info...and keep only a subset of vars
transition_points <- transition_points %>%
left_join(flt_rt, by = "flight_id") %>%
select(-c(ssr_codes, period_start, period_finish, source_ids, n))
write_csv(transition_points,
str_glue("data/egll_transition_points_{d}.csv", d = date_of_interest))
# join positions with relevant flight info...and keep only a subset of vars
pos_rt <- pos_rt_40NM %>%
left_join(flt_rt, by = "flight_id") %>%
select(flight_id,
timestamp,
latitude,
longitude,
altitude,
speed_gnd,
track_gnd,
vert_speed,
# on_ground,
# distance,
# distance_arp,
callsign,
aircraft_reg,
aircraft_type,
aircraft_address,
adep,
ades)
### WORK IN PROGRESS: START ###
pos_rt_40NM_sf <- pos_rt %>%
dplyr::group_by(flight_id) %>%
dplyr::arrange(timestamp) %>%
dplyr::filter(row_number() == 1) %>%
add_bearing(arp) %>%
add_quadrant() %>%
st_as_sf(
coords = c("longitude", "latitude"),
dim = "XY",
crs = 4326) %>%
ungroup() %>%
select(flight_id, timestamp, altitude, bearing, quadrant)
# TODO: see "Programming with dplyr"
# https://dplyr.tidyverse.org/articles/programming.html
# plot to check
add_vors <- function(egll_vor) {
ggplot2::geom_point(data = egll_vor,
ggplot2::aes(x = longitude, y = latitude),
colour = "blue",
size = 2)
}
add_label_vors <- function(egll_vor) {
displacement <- tribble(
~id, ~hjust, ~vjust,
"BNN", -0.2, 1.5,
"BIG", 1.2, -0.7,
"LAM", 1.2, 1.5,
"OCK", -0.2, -0.7
)
egll_vor <- egll_vor %>% left_join(displacement)
ggplot2::geom_text(data = egll_vor,
ggplot2::aes(x = longitude,
y = latitude,
label = id,
hjust = hjust,
vjust = vjust))
}
add_runways <- function(rw_sf) {
ggplot2::geom_sf(data = rw_sf,
size = 1.2,
alpha = 0)
}
add_holding_tracks <- function(holdings_sf, colour = "black") {
ggplot2::geom_sf(
data = holdings_sf,
colour = colour,
alpha = 0)
}
### WORK IN PROGRESS: END ###
f_ids <- pos_rt %>%
dplyr::distinct(flight_id) %>%
pull()
pos_rt.sf <- pos_rt %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
# In order to intercept trajectory portions within holding polygons
# we have to box altitudes. This is expecially critical for BIG
# AIP for BIG, LAM, OCK, BNN states minumum level at 7000 ft,
# so we should consider only reports somewhat above that, say 6700 ft
tracks <- pos_rt.sf %>%
filter(altitude > 6700) %>%
dplyr::group_by(flight_id) %>%
dplyr::arrange(timestamp) %>%
dplyr::summarise(geometry = st_combine(geometry)) %>%
sf::st_cast("LINESTRING") %>%
ungroup() %>%
mutate(flown_40NM = st_length(.))
# NOTE: there are strange trajectories that continue beyond the runway:
# these seems to be due to faulty CPRs
tma_times <- pos_rt.sf %>%
dplyr::group_by(flight_id) %>%
dplyr::arrange(timestamp) %>%
mutate(tma_entry_timestamp = first(timestamp), tma_exit_timestamp = last(timestamp),
time_tma = tma_exit_timestamp - tma_entry_timestamp) %>%
filter(row_number()==1) %>%
select(flight_id, time_tma, tma_entry_timestamp, tma_exit_timestamp) %>%
st_set_geometry(NULL)
details <- pos_rt.sf %>%
group_by(flight_id) %>%
filter(row_number() == 1) %>%
ungroup() %>%
st_set_geometry(NULL) %>%
select(-altitude, -speed_gnd, -track_gnd, -vert_speed)
# tracks <- tracks %>%
# left_join(details)
# ggplot2::ggplot() +
# add_runways(egll_rw_sf) +
# add_vors(egll_vor) +
# add_label_vors(egll_vor) +
# add_holding_tracks(hss_sf) +
# add_holding_tracks(holdings, colour = "red")
holding_tracks <- tracks %>%
st_intersection(holdings)
holding_distances <- holding_tracks %>%
mutate(flown_holding = st_length(.)) %>%
st_set_geometry(NULL)
holding_times <- pos_rt.sf %>%
filter(altitude > 6700) %>%
st_intersection(holdings) %>%
st_set_geometry(NULL) %>%
dplyr::group_by(flight_id) %>%
dplyr::arrange(timestamp) %>%
mutate(
holding_entry_timestamp = first(timestamp),
holding_exit_timestamp = last(timestamp),
time_holding = holding_exit_timestamp - holding_entry_timestamp) %>%
filter(row_number()==1) %>%
ungroup() %>%
select(flight_id, time_holding, holding_entry_timestamp, holding_exit_timestamp)
plot_tracks <- function(tks, from, to) {
ts <- tks %>%
filter(dplyr::row_number() >= from, dplyr::row_number() < to)
ggplot2::ggplot() +
add_runways(egll_rw_sf) +
add_vors(egll_vor) +
add_label_vors(egll_vor) +
geom_sf(data = ts, alpha = 1e-6) +
add_holding_tracks(hss_sf, colour = "yellow") +
add_holding_tracks(holdings, colour = "red")
}
# # visualize the different LINESTRINGs
# # (outcome from self-intersections)
# holding_tracks %>%
# filter(row_number() == 3) %>%
# st_cast("LINESTRING") %>%
# mutate(seq = row_number()) %>%
# ggplot() + geom_sf(aes(colour = seq))
# calc flown lengths within holding stacks
holding_tracks <- holding_tracks %>%
mutate(holding_length = st_length(.)) %>%
mutate(holding_length = drop_units(holding_length))
# # facet distribution of flown lengths by holding stack
holding_tracks %>%
ggplot(aes(x = holding_length, group = id, colour = id)) +
geom_density() +
facet_wrap(vars(id))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.