eneller/opensky-statistics
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

♻️ refactor all logic to main.rmd

Browse Source
This commit is contained in:
lukasadrion
2026-01-20 16:53:19 +01:00
parent eb49746268
commit 696f52eda3
2 changed files with 224 additions and 547 deletions
Download Patch File Download Diff File Expand all files Collapse all files

305
src/app.Rmd
View File

@@ -1,22 +1,13 @@
```{r backend, child="./main.Rmd"} ```{r backend, include=FALSE}
source("./main.Rmd") # Load all functions from main.Rmd
exists("getRouteSummary") knitr::purl("./main.Rmd", output = tempfile(), quiet = TRUE) |> source()
``` ```
# Web Interface # Web Interface
```{r shiny} ```{r shiny}
# Flight Trajectory Analysis - Shiny GUI Application # Flight Trajectory Analysis - Shiny GUI Application
# This app allows interactive selection of flights and displays trajectory analysis # This app allows interactive selection of flights and displays trajectory analysis
# All core functions are loaded from main.Rmd
library(shiny)
library(dplyr)
library(lubridate)
library(openSkies)
library(dotenv)
library(httr)
library(jsonlite)
library(trajr)
# UI Definition # UI Definition
ui <- fluidPage( ui <- fluidPage(
@@ -132,6 +123,7 @@ server <- function(input, output, session) {
status("Loading departures...") status("Loading departures...")
tryCatch({ tryCatch({
# Use getCredentials from main.Rmd
rv$creds <- getCredentials( rv$creds <- getCredentials(
client_id = input$client_id, client_id = input$client_id,
client_secret = input$client_secret client_secret = input$client_secret
@@ -199,28 +191,18 @@ server <- function(input, output, session) {
rv$current_icao <- icao24 rv$current_icao <- icao24
# Get track data # Use getAircraftTrack from main.Rmd
query <- list(icao24 = icao24, time = as.numeric(dep_time)) route_df <- getAircraftTrack(icao24, dep_time, rv$creds)
response <- makeAuthenticatedRequest('tracks/all', query, rv$creds)
if (httr::status_code(response) != 200) { if (is.null(route_df) || nrow(route_df) < 2) {
status(paste("Track data not available for", icao24, "(HTTP", httr::status_code(response), ")"))
return()
}
track_data <- fromJSON(content(response, as = "text", encoding = "UTF-8"))
if (is.null(track_data$path) || length(track_data$path) < 2) {
status(paste("No path data available for", icao24)) status(paste("No path data available for", icao24))
return() return()
} }
route_df <- as.data.frame(track_data$path)
colnames(route_df) <- c("time", "lat", "lon", "alt", "heading", "on_ground")
rv$current_route <- route_df rv$current_route <- route_df
# Create trajectory object # Use getTrajFromRoute from main.Rmd
rv$current_trj <- createTrajFromRoute(route_df) rv$current_trj <- getTrajFromRoute(route_df)
status(paste("Successfully analyzed", icao24, "with", nrow(route_df), "points")) status(paste("Successfully analyzed", icao24, "with", nrow(route_df), "points"))
# Switch to analysis tab # Switch to analysis tab
@@ -256,15 +238,10 @@ server <- function(input, output, session) {
# Characteristics table # Characteristics table
output$characteristics_table <- renderTable({ output$characteristics_table <- renderTable({
req(rv$current_trj) req(rv$current_trj)
calculateTrajectoryStats(rv$current_trj, format = "table")
trj <- rv$current_trj
data.frame <- calculateRouteCharacteristics(trj)
data.frame
}) })
# Batch analysis # Batch analysis
# FIXME use multiple flights from one aircraft instead of random flights of random aircrafts
observeEvent(input$batch_analyze, { observeEvent(input$batch_analyze, {
req(rv$departures, rv$creds) req(rv$departures, rv$creds)
@@ -284,47 +261,8 @@ server <- function(input, output, session) {
if (is.null(dep_time)) next if (is.null(dep_time)) next
params <- tryCatch({ # Use calculate_trajectory_params from main.Rmd
query <- list(icao24 = icao24, time = as.numeric(dep_time)) params <- calculate_trajectory_params(icao24, dep_time, rv$creds)
response <- makeAuthenticatedRequest('tracks/all', query, rv$creds)
if (httr::status_code(response) != 200) return(NULL)
track_data <- fromJSON(content(response, as = "text", encoding = "UTF-8"))
if (is.null(track_data$path) || length(track_data$path) < 3) return(NULL)
route_df <- as.data.frame(track_data$path)
colnames(route_df) <- c("time", "lat", "lon", "alt", "heading", "on_ground")
trj <- createTrajFromRoute(route_df)
duration <- TrajDuration(trj)
path_length <- TrajLength(trj)
diffusion_dist <- TrajDistance(trj)
straight <- TrajStraightness(trj)
mean_vel <- path_length / duration
fractal <- tryCatch({
min_step <- path_length / 100
max_step <- path_length / 2
if (min_step > 0 && max_step > min_step) {
step_sizes <- exp(seq(log(min_step), log(max_step), length.out = 10))
TrajFractalDimension(trj, stepSizes = step_sizes)
} else {
NA
}
}, error = function(e) NA)
data.frame(
icao24 = icao24,
diffusion_distance_km = diffusion_dist / 1000,
straightness = straight,
duration_min = duration / 60,
mean_velocity_kmh = mean_vel * 3.6,
fractal_dimension = fractal
)
}, error = function(e) NULL)
if (!is.null(params)) { if (!is.null(params)) {
all_trajectories[[length(all_trajectories) + 1]] <- params all_trajectories[[length(all_trajectories) + 1]] <- params
@@ -347,242 +285,37 @@ server <- function(input, output, session) {
}) })
}) })
# Statistics summary table # Statistics summary table - use calculateStatsSummary from main.Rmd
output$stats_summary_table <- renderTable({ output$stats_summary_table <- renderTable({
req(rv$trajectory_stats_df) req(rv$trajectory_stats_df)
calculateStatsSummary(rv$trajectory_stats_df) calculateStatsSummary(rv$trajectory_stats_df)
}) })
# Boxplots # Boxplots - use createBoxplots from main.Rmd
output$boxplots <- renderPlot({ output$boxplots <- renderPlot({
req(rv$trajectory_stats_df) req(rv$trajectory_stats_df)
createBoxplots(rv$trajectory_stats_df) createBoxplots(rv$trajectory_stats_df)
}) })
# Density plots # Density plots - use createDensityPlots from main.Rmd
output$density_plots <- renderPlot({ output$density_plots <- renderPlot({
req(rv$trajectory_stats_df) req(rv$trajectory_stats_df)
createDensityPlots(rv$trajectory_stats_df) createDensityPlots(rv$trajectory_stats_df)
}) })
# Histograms # Histograms - use createHistograms from main.Rmd
output$histograms <- renderPlot({ output$histograms <- renderPlot({
req(rv$trajectory_stats_df) req(rv$trajectory_stats_df)
createHistograms(rv$trajectory_stats_df) createHistograms(rv$trajectory_stats_df)
}) })
# Interpretation text # Interpretation text - use generateInterpretation from main.Rmd
output$interpretation_text <- renderText({ output$interpretation_text <- renderText({
req(rv$trajectory_stats_df) req(rv$trajectory_stats_df)
generateInterpretation(rv$trajectory_stats_df) generateInterpretation(rv$trajectory_stats_df)
}) })
} }
# Helper function to get parameter names and labels
getTrajectoryParams <- function() {
list(
params = c("diffusion_distance_km", "straightness", "duration_min",
"mean_velocity_kmh", "fractal_dimension"),
labels = c("Diffusion Distance (km)", "Straightness", "Duration (min)",
"Mean Velocity (km/h)", "Fractal Dimension")
)
}
# Calculate statistics summary table
calculateStatsSummary <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
stats_list <- lapply(seq_along(p$params), function(i) {
x <- trajectory_stats_df[[p$params[i]]]
x <- x[!is.na(x)]
if (length(x) < 2) return(NULL)
data.frame(
Parameter = p$labels[i],
N = length(x),
Mean = round(mean(x), 4),
Variance = round(var(x), 4),
Std_Dev = round(sd(x), 4),
Q1 = round(quantile(x, 0.25), 4),
Median = round(median(x), 4),
Q3 = round(quantile(x, 0.75), 4)
)
})
do.call(rbind, stats_list[!sapply(stats_list, is.null)])
}
# Create boxplots for trajectory statistics
createBoxplots <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
par(mfrow = c(2, 3))
for (i in seq_along(p$params)) {
data <- trajectory_stats_df[[p$params[i]]][!is.na(trajectory_stats_df[[p$params[i]]])]
if (length(data) >= 2) {
boxplot(data, main = p$labels[i], ylab = p$labels[i], col = "lightblue", border = "darkblue")
points(1, mean(data), pch = 18, col = "red", cex = 1.5)
}
}
par(mfrow = c(1, 1))
}
# Create density plots for trajectory statistics
createDensityPlots <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
par(mfrow = c(2, 3))
for (i in seq_along(p$params)) {
data <- trajectory_stats_df[[p$params[i]]][!is.na(trajectory_stats_df[[p$params[i]]])]
if (length(data) >= 3) {
dens <- density(data)
plot(dens, main = paste("Density:", p$labels[i]), xlab = p$labels[i], col = "darkblue", lwd = 2)
polygon(dens, col = rgb(0, 0, 1, 0.3), border = "darkblue")
abline(v = mean(data), col = "red", lwd = 2, lty = 2)
abline(v = median(data), col = "green", lwd = 2, lty = 3)
}
}
par(mfrow = c(1, 1))
}
# Create histograms for trajectory statistics
createHistograms <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
par(mfrow = c(2, 3))
for (i in seq_along(p$params)) {
data <- trajectory_stats_df[[p$params[i]]][!is.na(trajectory_stats_df[[p$params[i]]])]
if (length(data) >= 3) {
hist(data, probability = TRUE, main = paste("Histogram:", p$labels[i]),
xlab = p$labels[i], col = "lightgray", border = "darkgray")
lines(density(data), col = "red", lwd = 2)
}
}
par(mfrow = c(1, 1))
}
# Generate interpretation text for trajectory statistics
generateInterpretation <- function(trajectory_stats_df) {
df <- trajectory_stats_df
text <- "========== INTERPRETATION OF TRAJECTORY PARAMETERS ==========\n\n"
# Diffusion Distance
dd <- df$diffusion_distance_km[!is.na(df$diffusion_distance_km)]
if (length(dd) >= 2) {
text <- paste0(text, "1. DIFFUSION DISTANCE (Net Displacement):\n")
text <- paste0(text, " - Mean: ", round(mean(dd), 2), " km\n")
text <- paste0(text, " - Represents straight-line distance from origin to destination.\n")
text <- paste0(text, " - Variance: ", round(var(dd), 2), " (indicates diversity in flight distances)\n\n")
}
# Straightness
st <- df$straightness[!is.na(df$straightness)]
if (length(st) >= 2) {
text <- paste0(text, "2. STRAIGHTNESS INDEX:\n")
text <- paste0(text, " - Mean: ", round(mean(st), 4), " (range 0-1, where 1 = perfectly straight)\n")
text <- paste0(text, " - Values close to 1 indicate efficient, direct flight paths.\n")
text <- paste0(text, " - Lower values suggest deviations due to weather, airspace, or routing.\n\n")
}
# Duration
dur <- df$duration_min[!is.na(df$duration_min)]
if (length(dur) >= 2) {
text <- paste0(text, "3. DURATION OF TRAVEL:\n")
text <- paste0(text, " - Mean: ", round(mean(dur), 2), " minutes\n")
text <- paste0(text, " - Range: ", round(min(dur), 2), " - ", round(max(dur), 2), " minutes\n")
text <- paste0(text, " - IQR: ", round(IQR(dur), 2), " minutes (middle 50% of flights)\n\n")
}
# Velocity
vel <- df$mean_velocity_kmh[!is.na(df$mean_velocity_kmh)]
if (length(vel) >= 2) {
text <- paste0(text, "4. MEAN TRAVEL VELOCITY:\n")
text <- paste0(text, " - Mean: ", round(mean(vel), 2), " km/h\n")
text <- paste0(text, " - Typical commercial aircraft cruise: 800-900 km/h\n")
text <- paste0(text, " - Lower values may include taxi, takeoff, and landing phases.\n\n")
}
# Fractal Dimension
fd <- df$fractal_dimension[!is.na(df$fractal_dimension)]
if (length(fd) >= 2) {
text <- paste0(text, "5. FRACTAL DIMENSION:\n")
text <- paste0(text, " - Mean: ", round(mean(fd), 4), "\n")
text <- paste0(text, " - Value of 1.0 = perfectly straight line\n")
text <- paste0(text, " - Values closer to 2.0 = more complex, space-filling paths\n")
text <- paste0(text, " - Aircraft typically show low fractal dimension (efficient paths).\n\n")
}
text <- paste0(text, "========== END OF ANALYSIS ==========")
text
}
createTrajFromRoute <- function(route_df) {
tryCatch({
lat_ref <- route_df$lat[1]
lon_ref <- route_df$lon[1]
meters_per_deg_lat <- 111320
meters_per_deg_lon <- 111320 * cos(lat_ref * pi / 180)
x_meters <- (route_df$lon - lon_ref) * meters_per_deg_lon
y_meters <- (route_df$lat - lat_ref) * meters_per_deg_lat
time_seconds <- route_df$time - route_df$time[1]
trj <- TrajFromCoords(
data.frame(x = x_meters, y = y_meters, time = time_seconds),
xCol = "x", yCol = "y", timeCol = "time"
)
return(trj)
}, error = function(e) {
status(paste("Error creating trajectory object:", e$message))
})
}
calculateRouteCharacteristics <- function(trj) {
duration <- TrajDuration(trj)
path_length <- TrajLength(trj)
diffusion_distance <- TrajDistance(trj)
straightness <- TrajStraightness(trj)
mean_velocity <- path_length / duration
fractal_dim <- tryCatch({
min_step <- path_length / 100
max_step <- path_length / 2
if (min_step > 0 && max_step > min_step) {
step_sizes <- exp(seq(log(min_step), log(max_step), length.out = 10))
TrajFractalDimension(trj, stepSizes = step_sizes)
} else {
NA
}
}, error = function(e) NA)
return (data.frame(
Parameter = c(
"Duration (s)", "Duration (min)",
"Path Length (km)",
"Duffusion Distance (m)",
"Diffusion Distance (km)",
"Straightness Index",
"Mean Velocity (km/h)",
"Fractal Dimension"
),
Value = c(
duration_s = round(duration, 2),
duration_min = round(duration / 60, 2),
path_length_km = round(path_length / 1000, 2),
diffusion_distance_m = round(diffusion_distance, 2),
diffusion_distance_km = round(diffusion_distance / 1000, 2),
straightness_index = round(straightness, 4),
mean_velocity_kmh = round(mean_velocity *3.6, 2),
fractal_dimension = round(fractal_dim, 4)
)
)
)
}
# Run the application # Run the application
shinyApp(ui = ui, server = server) shinyApp(ui = ui, server = server)
``` ```

466
src/main.Rmd
View File

@@ -1,5 +1,5 @@
--- ---
title: "Topic 8" title: "Topic 8 - Flight Trajectory Analysis"
output: output:
pdf_document: default pdf_document: default
html_document: default html_document: default
@@ -21,6 +21,7 @@ library(dotenv)
library(httr) library(httr)
library(jsonlite) library(jsonlite)
library(trajr) library(trajr)
library(shiny)
``` ```
# Download flights # Download flights
@@ -38,373 +39,316 @@ getFlights <- function(icao, time, creds){
flights <-getAircraftFlights(icao, startTime = time - days(1), endTime = time, credentials = creds ) flights <-getAircraftFlights(icao, startTime = time - days(1), endTime = time, credentials = creds )
return(flights) return(flights)
} }
icao <- departures[[1]][["ICAO24"]]
flights <- getFlights(icao,Sys.time(), creds)
# TODO map from all available flights to tracks
query <- list(icao24= icao, time=as.numeric(flights[[1]][["departure_time"]]))
# make a manual request because this API endpoint is considered experimental # Get aircraft track from OpenSky API
getAircraftTrack <- function(icao, time, creds){ getAircraftTrack <- function(icao, time, creds) {
query <- list(icao24= icao, time=as.numeric(time)) query <- list(icao24 = icao, time = as.numeric(time))
response <-makeAuthenticatedRequest('tracks/all',query, creds) response <- makeAuthenticatedRequest('tracks/all', query, creds)
track_data <- fromJSON(content(response, as = "text", encoding = "UTF-8")) track_data <- fromJSON(content(response, as = "text", encoding = "UTF-8"))
if (!is.null(track_data$path) && length(track_data$path) > 0) { if (!is.null(track_data$path) && length(track_data$path) > 0) {
route_df <- as.data.frame(track_data$path) route_df <- as.data.frame(track_data$path)
colnames(route_df) <- c("time", "lat", "lon", "alt", "heading", "on_ground") colnames(route_df) <- c("time", "lat", "lon", "alt", "heading", "on_ground")
return(route_df) return(route_df)
} }
return(NULL) return(NULL)
} }
route_df <- getAircraftTrack(icao, time=flights[[1]][["departure_time"]], creds = creds)
if(!is.null(route_df)){
message("Loading of route successfull! Points: ", nrow(route_df))
plot(route_df$lon, route_df$lat, type="o", pch=20, col="blue",
main=paste("Geographic route of", icao),
xlab="Longitude", ylab="Latitude")
plot(route_df$time, route_df$alt, type="l", col="red", lwd=2,
main=paste("Altitude profile of", icao),
xlab="Time (Unix)", ylab="Height (Meter)")
} else {
print("No path points from api")
}
``` ```
## Trajectory Conversion Functions
# Trajectory Characteristics Analysis ```{r trajectory-functions}
# Convert route to distance in meters
## Distance Approximation getRouteDistance <- function(route_df) {
```{r traj-dist}
getRouteDistance<- function(route_df){
# Convert lat/lon to approximate meters (using simple equirectangular projection)
# Reference point: first coordinate
lat_ref <- route_df$lat[1] lat_ref <- route_df$lat[1]
lon_ref <- route_df$lon[1] lon_ref <- route_df$lon[1]
# Convert to meters (approximate)
meters_per_deg_lat <- 111320 meters_per_deg_lat <- 111320
meters_per_deg_lon <- 111320 * cos(lat_ref * pi / 180) meters_per_deg_lon <- 111320 * cos(lat_ref * pi / 180)
x_meters <- (route_df$lon - lon_ref) * meters_per_deg_lon x_meters <- (route_df$lon - lon_ref) * meters_per_deg_lon
y_meters <- (route_df$lat - lat_ref) * meters_per_deg_lat y_meters <- (route_df$lat - lat_ref) * meters_per_deg_lat
return(list('x' = x_meters, 'y' = y_meters)) return(list('x' = x_meters, 'y' = y_meters))
} }
getRouteTime <- function(route_df){ # Get time in seconds from start
getRouteTime <- function(route_df) {
return(route_df$time - route_df$time[1]) return(route_df$time - route_df$time[1])
} }
getTrajFromRoute <- function(route_df){ # Create trajr object from route
getTrajFromRoute <- function(route_df) {
meters <- getRouteDistance(route_df) meters <- getRouteDistance(route_df)
time <- getRouteTime(route_df) time <- getRouteTime(route_df)
trj <- TrajFromCoords( trj <- TrajFromCoords(
data.frame(x = meters$x, y = meters$y, time = time), data.frame(x = meters$x, y = meters$y, time = time),
xCol = "x", yCol = "y", timeCol = "time" xCol = "x", yCol = "y", timeCol = "time"
) )
return(trj)
} }
getRouteSummary<-function(route_df, icao){
meters <- getRouteDistance(route_df) # Calculate trajectory characteristics
x_meters <- meters$x # Input: either route_df (data.frame with lat/lon) or trj (trajr object)
y_meters <- meters$y # format: "row" for batch analysis (one row per flight), "table" for single flight display
time_seconds <- getRouteTime(route_df) # FIXME for batch analysis: use the same aircraft
calculateTrajectoryStats <- function(input, icao = NULL, format = "row") {
# Create trajr trajectory object # Determine if input is route_df or trj
trj <- getTrajFromRoute(route_df) if (inherits(input, "Trajectory")) {
trj <- input
# Calculate trajectory characteristics } else {
trj <- getTrajFromRoute(input)
# 1. Duration of travel (seconds) }
duration <- TrajDuration(trj)
# Calculate all metrics
# 2. Total path length (meters) duration <- TrajDuration(trj)
path_length <- TrajLength(trj) path_length <- TrajLength(trj)
diffusion_distance <- TrajDistance(trj)
# 3. Diffusion distance (net displacement - straight line from start to end) straightness <- TrajStraightness(trj)
diffusion_distance <- TrajDistance(trj) mean_velocity <- path_length / duration
# 4. Straightness index (ratio of net displacement to path length, 0-1) fractal_dim <- tryCatch({
straightness <- TrajStraightness(trj) min_step <- path_length / 100
max_step <- path_length / 2
# 5. Mean travel velocity (meters/second) if (min_step > 0 && max_step > min_step) {
mean_velocity <- path_length / duration
# 6. Fractal dimension (using divider method)
# Note: requires sufficient points for accurate estimation
fractal_dim <- tryCatch({
# Calculate appropriate step sizes based on trajectory length
min_step <- TrajLength(trj) / 100
max_step <- TrajLength(trj) / 2
step_sizes <- exp(seq(log(min_step), log(max_step), length.out = 10)) step_sizes <- exp(seq(log(min_step), log(max_step), length.out = 10))
TrajFractalDimension(trj, stepSizes = step_sizes) TrajFractalDimension(trj, stepSizes = step_sizes)
}, error = function(e) { } else {
message("Could not calculate fractal dimension: ", e$message)
NA NA
}) }
}, error = function(e) NA)
# Return format based on use case
if (format == "table") {
# For single flight display (Parameter | Value)
return(data.frame(
Parameter = c(
"Duration (s)", "Duration (min)",
"Path Length (km)",
"Diffusion Distance (m)",
"Diffusion Distance (km)",
"Straightness Index",
"Mean Velocity (km/h)",
"Fractal Dimension"
),
Value = c(
round(duration, 2),
round(duration / 60, 2),
round(path_length / 1000, 2),
round(diffusion_distance, 2),
round(diffusion_distance / 1000, 2),
round(straightness, 4),
round(mean_velocity * 3.6, 2),
round(fractal_dim, 4)
)
))
} else {
# For batch analysis (one row per flight)
return(data.frame( return(data.frame(
icao24 = icao, icao24 = icao,
diffusion_distance_km = diffusion_distance / 1000, # meters to kilometers diffusion_distance_km = diffusion_distance / 1000,
path_length_km = path_length / 1000, # meters to kilometers path_length_km = path_length / 1000,
straightness = straightness, straightness = straightness,
duration_mins = duration / 60, # seconds to minutes duration_min = duration / 60,
mean_velocity_kmh = mean_velocity * 3.6, mean_velocity_kmh = mean_velocity * 3.6,
fractal_dimension = fractal_dim fractal_dimension = fractal_dim
)) ))
}
} }
print(getRouteSummary(route_df, icao)) # Calculate trajectory parameters for a single flight
trj <- getTrajFromRoute(route_df)
plot(trj, main = paste("Trajectory of", icao))
```
# Statistical Analysis of Multiple Trajectories
```{r multi-trajectory-analysis}
# Function to calculate trajectory characteristics for a single flight
calculate_trajectory_params <- function(icao, departure_time, creds) { calculate_trajectory_params <- function(icao, departure_time, creds) {
tryCatch({ tryCatch({
route_df <- getAircraftTrack(icao, departure_time, creds) route_df <- getAircraftTrack(icao, departure_time, creds)
if (is.null(route_df) || nrow(route_df) < 3) return(NULL) if (is.null(route_df) || nrow(route_df) < 3) return(NULL)
# Fractal dimension return(calculateTrajectoryStats(route_df, icao = icao, format = "row"))
fractal <- tryCatch({
trj <- getTrajFromRoute(route_df)
path_length <- TrajLength(trj)
min_step <- path_length / 100
max_step <- path_length / 2
if (min_step > 0 && max_step > min_step) {
step_sizes <- exp(seq(log(min_step), log(max_step), length.out = 10))
TrajFractalDimension(trj, stepSizes = step_sizes)
} else {
NA
}
}, error = function(e) NA)
return(getRouteSummary(route_df, icao))
}, error = function(e) { }, error = function(e) {
message("Error processing ", icao, ": ", e$message) message("Error processing ", icao, ": ", e$message)
return(NULL) return(NULL)
}) })
} }
# Collect trajectory data from multiple departures
message("Collecting trajectory data from departures...")
all_trajectories <- list()
# Process available departures (limit to avoid API rate limits)
n_departures <- min(length(departures), 20)
for (i in 1:n_departures) {
dep <- departures[[i]]
icao24 <- dep[["ICAO24"]]
dep_time <- dep[["departure_time"]] # Use departure time directly from departures list
# Skip if no departure time available
if (is.null(dep_time)) {
message("Skipping ", icao24, ": no departure time")
next
}
params <- calculate_trajectory_params(icao24, dep_time, creds)
if (!is.null(params)) {
all_trajectories[[length(all_trajectories) + 1]] <- params
message("Successfully processed trajectory for ", icao24)
}
Sys.sleep(0.5) # Rate limiting
}
# Combine all trajectory data
if (length(all_trajectories) > 0) {
trajectory_stats_df <- do.call(rbind, all_trajectories)
message("Successfully collected ", nrow(trajectory_stats_df), " trajectories")
print(trajectory_stats_df)
} else {
message("No trajectory data collected")
}
``` ```
# Basic Statistical Analysis of Trajectory Parameters ## Statistical Helper Functions
```{r statistical-analysis} ```{r stat-functions}
if (exists("trajectory_stats_df") && nrow(trajectory_stats_df) >= 2) { # Get parameter names and labels for trajectory statistics
getTrajectoryParams <- function() {
list(
params = c("diffusion_distance_km", "straightness", "duration_min",
"mean_velocity_kmh", "fractal_dimension"),
labels = c("Diffusion Distance (km)", "Straightness", "Duration (min)",
"Mean Velocity (km/h)", "Fractal Dimension")
)
}
# Calculate statistics summary table
calculateStatsSummary <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
# Parameters to analyze stats_list <- lapply(seq_along(p$params), function(i) {
params_to_analyze <- c("diffusion_distance_km", "straightness", "duration_mins", x <- trajectory_stats_df[[p$params[i]]]
"mean_velocity_kmh", "fractal_dimension")
param_labels <- c("Diffusion Distance (km)", "Straightness Index",
"Duration (min)", "Mean Velocity (km/h)", "Fractal Dimension")
# Function to calculate comprehensive statistics
calc_stats <- function(x, param_name) {
x <- x[!is.na(x)] x <- x[!is.na(x)]
if (length(x) < 2) return(NULL) if (length(x) < 2) return(NULL)
data.frame( data.frame(
Parameter = param_name, Parameter = p$labels[i],
N = length(x), N = length(x),
Mean = round(mean(x), 4), Mean = round(mean(x), 4),
Variance = round(var(x), 4), Variance = round(var(x), 4),
Std_Dev = round(sd(x), 4), Std_Dev = round(sd(x), 4),
Min = round(min(x), 4),
Q1 = round(quantile(x, 0.25), 4), Q1 = round(quantile(x, 0.25), 4),
Median = round(median(x), 4), Median = round(median(x), 4),
Q3 = round(quantile(x, 0.75), 4), Q3 = round(quantile(x, 0.75), 4)
Max = round(max(x), 4),
IQR = round(IQR(x), 4)
) )
} })
# Calculate statistics for each parameter do.call(rbind, stats_list[!sapply(stats_list, is.null)])
stats_list <- list() }
for (i in seq_along(params_to_analyze)) { ```
param <- params_to_analyze[i]
label <- param_labels[i] ## Visualization Functions
stats_list[[i]] <- calc_stats(trajectory_stats_df[[param]], label) ```{r viz-functions}
} # Create boxplots for trajectory statistics
createBoxplots <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
# Combine into summary table
stats_summary <- do.call(rbind, stats_list[!sapply(stats_list, is.null)])
cat("\n========== STATISTICAL SUMMARY ==========\n\n")
print(stats_summary, row.names = FALSE)
# Boxplots for each parameter
par(mfrow = c(2, 3)) par(mfrow = c(2, 3))
for (i in seq_along(p$params)) {
for (i in seq_along(params_to_analyze)) { data <- trajectory_stats_df[[p$params[i]]][!is.na(trajectory_stats_df[[p$params[i]]])]
param <- params_to_analyze[i]
label <- param_labels[i]
data <- trajectory_stats_df[[param]][!is.na(trajectory_stats_df[[param]])]
if (length(data) >= 2) { if (length(data) >= 2) {
boxplot(data, boxplot(data, main = p$labels[i], ylab = p$labels[i], col = "lightblue", border = "darkblue")
main = paste("Boxplot:", label),
ylab = label,
col = "lightblue",
border = "darkblue")
# Add mean point
points(1, mean(data), pch = 18, col = "red", cex = 1.5) points(1, mean(data), pch = 18, col = "red", cex = 1.5)
} }
} }
par(mfrow = c(1, 1)) par(mfrow = c(1, 1))
}
# Create density plots for trajectory statistics
createDensityPlots <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
# Density plots for each parameter
par(mfrow = c(2, 3)) par(mfrow = c(2, 3))
for (i in seq_along(p$params)) {
for (i in seq_along(params_to_analyze)) { data <- trajectory_stats_df[[p$params[i]]][!is.na(trajectory_stats_df[[p$params[i]]])]
param <- params_to_analyze[i]
label <- param_labels[i]
data <- trajectory_stats_df[[param]][!is.na(trajectory_stats_df[[param]])]
if (length(data) >= 3) { if (length(data) >= 3) {
dens <- density(data, na.rm = TRUE) dens <- density(data)
plot(dens, plot(dens, main = paste("Density:", p$labels[i]), xlab = p$labels[i], col = "darkblue", lwd = 2)
main = paste("Density:", label),
xlab = label,
ylab = "Density",
col = "darkblue",
lwd = 2)
polygon(dens, col = rgb(0, 0, 1, 0.3), border = "darkblue") polygon(dens, col = rgb(0, 0, 1, 0.3), border = "darkblue")
# Add vertical lines for mean and median
abline(v = mean(data), col = "red", lwd = 2, lty = 2) abline(v = mean(data), col = "red", lwd = 2, lty = 2)
abline(v = median(data), col = "green", lwd = 2, lty = 3) abline(v = median(data), col = "green", lwd = 2, lty = 3)
legend("topright", legend = c("Mean", "Median"),
col = c("red", "green"), lty = c(2, 3), lwd = 2, cex = 0.7)
} }
} }
par(mfrow = c(1, 1)) par(mfrow = c(1, 1))
}
# Create histograms for trajectory statistics
createHistograms <- function(trajectory_stats_df) {
p <- getTrajectoryParams()
# Histogram with density overlay
par(mfrow = c(2, 3)) par(mfrow = c(2, 3))
for (i in seq_along(p$params)) {
for (i in seq_along(params_to_analyze)) { data <- trajectory_stats_df[[p$params[i]]][!is.na(trajectory_stats_df[[p$params[i]]])]
param <- params_to_analyze[i]
label <- param_labels[i]
data <- trajectory_stats_df[[param]][!is.na(trajectory_stats_df[[param]])]
if (length(data) >= 3) { if (length(data) >= 3) {
hist(data, hist(data, probability = TRUE, main = paste("Histogram:", p$labels[i]),
probability = TRUE, xlab = p$labels[i], col = "lightgray", border = "darkgray")
main = paste("Histogram:", label),
xlab = label,
col = "lightgray",
border = "darkgray")
# Overlay density curve
lines(density(data), col = "red", lwd = 2) lines(density(data), col = "red", lwd = 2)
} }
} }
par(mfrow = c(1, 1)) par(mfrow = c(1, 1))
} else {
message("Insufficient trajectory data for statistical analysis (need at least 2 trajectories)")
} }
```
# Interpretation of Results # Generate interpretation text for trajectory statistics
```{r interpretation} generateInterpretation <- function(trajectory_stats_df) {
if (exists("trajectory_stats_df") && nrow(trajectory_stats_df) >= 2) { df <- trajectory_stats_df
cat("\n========== INTERPRETATION OF TRAJECTORY PARAMETERS ==========\n\n") text <- "========== INTERPRETATION OF TRAJECTORY PARAMETERS ==========\n\n"
# Diffusion Distance dd <- df$diffusion_distance_km[!is.na(df$diffusion_distance_km)]
dd <- trajectory_stats_df$diffusion_distance_km[!is.na(trajectory_stats_df$diffusion_distance_km)]
if (length(dd) >= 2) { if (length(dd) >= 2) {
cat("1. DIFFUSION DISTANCE (Net Displacement):\n") text <- paste0(text, "1. DIFFUSION DISTANCE (Net Displacement):\n")
cat(" - Mean:", round(mean(dd), 2), "km\n") text <- paste0(text, " - Mean: ", round(mean(dd), 2), " km\n")
cat(" - This represents the straight-line distance from origin to destination.\n") text <- paste0(text, " - Represents straight-line distance from origin to destination.\n")
cat(" - High variance (", round(var(dd), 2), ") indicates diverse flight distances.\n\n") text <- paste0(text, " - Variance: ", round(var(dd), 2), " (indicates diversity in flight distances)\n\n")
} }
# Straightness st <- df$straightness[!is.na(df$straightness)]
st <- trajectory_stats_df$straightness[!is.na(trajectory_stats_df$straightness)]
if (length(st) >= 2) { if (length(st) >= 2) {
cat("2. STRAIGHTNESS INDEX:\n") text <- paste0(text, "2. STRAIGHTNESS INDEX:\n")
cat(" - Mean:", round(mean(st), 4), "(range 0-1, where 1 = perfectly straight)\n") text <- paste0(text, " - Mean: ", round(mean(st), 4), " (range 0-1, where 1 = perfectly straight)\n")
cat(" - Values close to 1 indicate efficient, direct flight paths.\n") text <- paste0(text, " - Values close to 1 indicate efficient, direct flight paths.\n")
cat(" - Lower values suggest deviations due to weather, airspace, or routing.\n\n") text <- paste0(text, " - Lower values suggest deviations due to weather, airspace, or routing.\n\n")
} }
# Duration dur <- df$duration_min[!is.na(df$duration_min)]
dur <- trajectory_stats_df$duration_min[!is.na(trajectory_stats_df$duration_min)]
if (length(dur) >= 2) { if (length(dur) >= 2) {
cat("3. DURATION OF TRAVEL:\n") text <- paste0(text, "3. DURATION OF TRAVEL:\n")
cat(" - Mean:", round(mean(dur), 2), "minutes\n") text <- paste0(text, " - Mean: ", round(mean(dur), 2), " minutes\n")
cat(" - Range:", round(min(dur), 2), "-", round(max(dur), 2), "minutes\n") text <- paste0(text, " - Range: ", round(min(dur), 2), " - ", round(max(dur), 2), " minutes\n")
cat(" - IQR:", round(IQR(dur), 2), "minutes (middle 50% of flights)\n\n") text <- paste0(text, " - IQR: ", round(IQR(dur), 2), " minutes (middle 50% of flights)\n\n")
} }
# Velocity vel <- df$mean_velocity_kmh[!is.na(df$mean_velocity_kmh)]
vel <- trajectory_stats_df$mean_velocity_kmh[!is.na(trajectory_stats_df$mean_velocity_kmh)]
if (length(vel) >= 2) { if (length(vel) >= 2) {
cat("4. MEAN TRAVEL VELOCITY:\n") text <- paste0(text, "4. MEAN TRAVEL VELOCITY:\n")
cat(" - Mean:", round(mean(vel), 2), "km/h\n") text <- paste0(text, " - Mean: ", round(mean(vel), 2), " km/h\n")
cat(" - Typical commercial aircraft cruise: 800-900 km/h\n") text <- paste0(text, " - Typical commercial aircraft cruise: 800-900 km/h\n")
cat(" - Lower values may include taxi, takeoff, and landing phases.\n\n") text <- paste0(text, " - Lower values may include taxi, takeoff, and landing phases.\n\n")
} }
# Fractal Dimension fd <- df$fractal_dimension[!is.na(df$fractal_dimension)]
fd <- trajectory_stats_df$fractal_dimension[!is.na(trajectory_stats_df$fractal_dimension)]
if (length(fd) >= 2) { if (length(fd) >= 2) {
cat("5. FRACTAL DIMENSION:\n") text <- paste0(text, "5. FRACTAL DIMENSION:\n")
cat(" - Mean:", round(mean(fd), 4), "\n") text <- paste0(text, " - Mean: ", round(mean(fd), 4), "\n")
cat(" - Value of 1.0 = perfectly straight line\n") text <- paste0(text, " - Value of 1.0 = perfectly straight line\n")
cat(" - Values closer to 2.0 = more complex, space-filling paths\n") text <- paste0(text, " - Values closer to 2.0 = more complex, space-filling paths\n")
cat(" - Aircraft typically show low fractal dimension (efficient paths).\n\n") text <- paste0(text, " - Aircraft typically show low fractal dimension (efficient paths).\n\n")
} }
cat("========== END OF ANALYSIS ==========\n") text <- paste0(text, "========== END OF ANALYSIS ==========")
text
}
```
# Example Usage (Demo)
```{r demo, eval=FALSE}
# This section shows how to use the functions above
# Set eval=TRUE to run this demo
# Get credentials
creds <- getCredentials()
# Get departures from Frankfurt airport
time_now <- Sys.time()
departures <- getAirportDepartures(
airport = "EDDF",
startTime = time_now - hours(1),
endTime = time_now,
credentials = creds
)
# Get first departure's track
if (length(departures) > 0) {
icao <- departures[[1]][["ICAO24"]]
dep_time <- departures[[1]][["departure_time"]]
route_df <- getAircraftTrack(icao, dep_time, creds)
if (!is.null(route_df)) {
# Plot route
plot(route_df$lon, route_df$lat, type = "o", pch = 20, col = "blue",
main = paste("Geographic route of", icao),
xlab = "Longitude", ylab = "Latitude")
# Plot altitude
plot(route_df$time, route_df$alt, type = "l", col = "red", lwd = 2,
main = paste("Altitude profile of", icao),
xlab = "Time (Unix)", ylab = "Height (Meter)")
# Get summary
print(getRouteSummary(route_df, icao))
# Plot trajectory
trj <- getTrajFromRoute(route_df)
plot(trj, main = paste("Trajectory of", icao))
}
} }
``` ```