diff --git a/src/app.R b/src/app.R new file mode 100644 index 0000000..e1da9b1 --- /dev/null +++ b/src/app.R @@ -0,0 +1,553 @@ +# Flight Trajectory Analysis - Shiny GUI Application +# This app allows interactive selection of flights and displays trajectory analysis + +library(shiny) +library(dplyr) +library(lubridate) +library(openSkies) +library(dotenv) +library(httr) +library(jsonlite) +library(trajr) + +# UI Definition +ui <- fluidPage( + titlePanel("Flight Trajectory Analysis - GUI"), + + sidebarLayout( + sidebarPanel( + width = 3, + + h4("OpenSky Credentials"), + textInput("client_id", "Client ID:", value = Sys.getenv('OPENSKY_CLIENT_ID')), + passwordInput("client_secret", "Client Secret:", value = Sys.getenv('OPENSKY_CLIENT_SECRET')), + + hr(), + + h4("Airport Selection"), + textInput("airport_code", "Airport ICAO Code:", value = "EDDF"), + sliderInput("hours_back", "Hours back from now:", min = 1, max = 12, value = 1), + + actionButton("load_departures", "Load Departures", class = "btn-primary"), + + hr(), + + h4("Flight Selection"), + selectInput("selected_flight", "Select Flight:", choices = NULL), + + actionButton("analyze_flight", "Analyze Selected Flight", class = "btn-success"), + + hr(), + + h4("Batch Analysis"), + numericInput("batch_size", "Number of flights to analyze:", value = 10, min = 2, max = 50), + actionButton("batch_analyze", "Run Batch Analysis", class = "btn-warning"), + + hr(), + + verbatimTextOutput("status_text") + ), + + mainPanel( + width = 9, + + tabsetPanel( + id = "main_tabs", + + tabPanel("Departures List", + h4("Available Departures"), + tableOutput("departures_table") + ), + + tabPanel("Single Flight Analysis", + fluidRow( + column(6, plotOutput("route_plot", height = "400px")), + column(6, plotOutput("altitude_plot", height = "400px")) + ), + fluidRow( + column(6, plotOutput("trajectory_plot", height = "400px")), + column(6, + h4("Trajectory Characteristics"), + tableOutput("characteristics_table")) + ) + ), + + tabPanel("Statistical Analysis", + h4("Multiple Trajectory Statistics"), + tableOutput("stats_summary_table"), + hr(), + fluidRow( + column(12, plotOutput("boxplots", height = "500px")) + ), + fluidRow( + column(12, plotOutput("density_plots", height = "500px")) + ), + fluidRow( + column(12, plotOutput("histograms", height = "500px")) + ) + ), + + tabPanel("Interpretation", + h4("Analysis Interpretation"), + verbatimTextOutput("interpretation_text") + ) + ) + ) + ) +) + +# Server Logic +server <- function(input, output, session) { + + # Reactive values to store data + rv <- reactiveValues( + creds = NULL, + departures = NULL, + departures_df = NULL, + current_route = NULL, + current_trj = NULL, + current_icao = NULL, + trajectory_stats_df = NULL + ) + + # Status message + status <- reactiveVal("Ready. Enter credentials and load departures.") + + output$status_text <- renderText({ + status() + }) + + # Load departures + observeEvent(input$load_departures, { + req(input$client_id, input$client_secret, input$airport_code) + + status("Loading departures...") + + tryCatch({ + rv$creds <- getCredentials( + client_id = input$client_id, + client_secret = input$client_secret + ) + + time_now <- Sys.time() + rv$departures <- getAirportDepartures( + airport = input$airport_code, + startTime = time_now - hours(input$hours_back), + endTime = time_now, + credentials = rv$creds + ) + + if (length(rv$departures) > 0) { + # Create departures dataframe for display + departures_list <- lapply(seq_along(rv$departures), function(i) { + dep <- rv$departures[[i]] + data.frame( + Index = i, + ICAO24 = dep[["ICAO24"]] %||% NA, + Callsign = dep[["callsign"]] %||% NA, + Origin = dep[["estDepartureAirport"]] %||% NA, + Destination = dep[["estArrivalAirport"]] %||% NA, + DepartureTime = as.POSIXct(dep[["departure_time"]] %||% NA, origin = "1970-01-01"), + stringsAsFactors = FALSE + ) + }) + rv$departures_df <- do.call(rbind, departures_list) + + # Update flight selection dropdown + choices <- setNames( + seq_along(rv$departures), + paste(rv$departures_df$ICAO24, "-", rv$departures_df$Callsign, + "(", rv$departures_df$Destination, ")") + ) + updateSelectInput(session, "selected_flight", choices = choices) + + status(paste("Loaded", length(rv$departures), "departures from", input$airport_code)) + } else { + status("No departures found for the selected time period.") + } + + }, error = function(e) { + status(paste("Error loading departures:", e$message)) + }) + }) + + # Display departures table + output$departures_table <- renderTable({ + req(rv$departures_df) + rv$departures_df + }) + + # Analyze selected flight + observeEvent(input$analyze_flight, { + req(rv$departures, input$selected_flight, rv$creds) + + status("Analyzing selected flight...") + + tryCatch({ + idx <- as.integer(input$selected_flight) + dep <- rv$departures[[idx]] + icao24 <- dep[["ICAO24"]] + dep_time <- dep[["departure_time"]] + + rv$current_icao <- icao24 + + # Get track data + query <- list(icao24 = icao24, time = as.numeric(dep_time)) + response <- makeAuthenticatedRequest('tracks/all', query, rv$creds) + + if (httr::status_code(response) != 200) { + 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)) + 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 + + # Create trajectory object + 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] + + rv$current_trj <- TrajFromCoords( + data.frame(x = x_meters, y = y_meters, time = time_seconds), + xCol = "x", yCol = "y", timeCol = "time" + ) + + status(paste("Successfully analyzed", icao24, "with", nrow(route_df), "points")) + + # Switch to analysis tab + updateTabsetPanel(session, "main_tabs", selected = "Single Flight Analysis") + + }, error = function(e) { + status(paste("Error analyzing flight:", e$message)) + }) + }) + + # Route plot + output$route_plot <- renderPlot({ + req(rv$current_route) + plot(rv$current_route$lon, rv$current_route$lat, type = "o", pch = 20, col = "blue", + main = paste("Geographic Route of", rv$current_icao), + xlab = "Longitude", ylab = "Latitude") + }) + + # Altitude plot + output$altitude_plot <- renderPlot({ + req(rv$current_route) + plot(rv$current_route$time, rv$current_route$alt, type = "l", col = "red", lwd = 2, + main = paste("Altitude Profile of", rv$current_icao), + xlab = "Time (Unix)", ylab = "Altitude (m)") + }) + + # Trajectory plot + output$trajectory_plot <- renderPlot({ + req(rv$current_trj) + plot(rv$current_trj, main = paste("Trajectory of", rv$current_icao)) + }) + + # Characteristics table + output$characteristics_table <- renderTable({ + req(rv$current_trj) + + trj <- rv$current_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) + + data.frame( + Parameter = c( + "Duration (s)", "Duration (min)", + "Path Length (m)", "Path Length (km)", + "Diffusion Distance (m)", "Diffusion Distance (km)", + "Straightness Index", + "Mean Velocity (m/s)", "Mean Velocity (km/h)", + "Fractal Dimension" + ), + Value = c( + round(duration, 2), round(duration / 60, 2), + round(path_length, 2), round(path_length / 1000, 2), + round(diffusion_distance, 2), round(diffusion_distance / 1000, 2), + round(straightness, 4), + round(mean_velocity, 2), round(mean_velocity * 3.6, 2), + round(fractal_dim, 4) + ) + ) + }) + + # Batch analysis + observeEvent(input$batch_analyze, { + req(rv$departures, rv$creds) + + status("Running batch analysis...") + + tryCatch({ + all_trajectories <- list() + n_departures <- min(length(rv$departures), input$batch_size) + + withProgress(message = 'Analyzing flights', value = 0, { + for (i in 1:n_departures) { + dep <- rv$departures[[i]] + icao24 <- dep[["ICAO24"]] + dep_time <- dep[["departure_time"]] + + incProgress(1/n_departures, detail = paste("Processing", icao24)) + + if (is.null(dep_time)) next + + params <- tryCatch({ + query <- list(icao24 = icao24, time = as.numeric(dep_time)) + 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") + + 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" + ) + + 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)) { + all_trajectories[[length(all_trajectories) + 1]] <- params + } + + Sys.sleep(0.3) + } + }) + + if (length(all_trajectories) > 0) { + rv$trajectory_stats_df <- do.call(rbind, all_trajectories) + status(paste("Batch analysis complete:", nrow(rv$trajectory_stats_df), "trajectories analyzed")) + updateTabsetPanel(session, "main_tabs", selected = "Statistical Analysis") + } else { + status("No trajectory data collected in batch analysis") + } + + }, error = function(e) { + status(paste("Error in batch analysis:", e$message)) + }) + }) + + # Statistics summary table + output$stats_summary_table <- renderTable({ + req(rv$trajectory_stats_df) + + 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") + + stats_list <- lapply(seq_along(params), function(i) { + x <- rv$trajectory_stats_df[[params[i]]] + x <- x[!is.na(x)] + if (length(x) < 2) return(NULL) + + data.frame( + Parameter = 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)]) + }) + + # Boxplots + output$boxplots <- renderPlot({ + req(rv$trajectory_stats_df) + + 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") + + par(mfrow = c(2, 3)) + for (i in seq_along(params)) { + data <- rv$trajectory_stats_df[[params[i]]][!is.na(rv$trajectory_stats_df[[params[i]]])] + if (length(data) >= 2) { + boxplot(data, main = labels[i], ylab = labels[i], col = "lightblue", border = "darkblue") + points(1, mean(data), pch = 18, col = "red", cex = 1.5) + } + } + par(mfrow = c(1, 1)) + }) + + # Density plots + output$density_plots <- renderPlot({ + req(rv$trajectory_stats_df) + + 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") + + par(mfrow = c(2, 3)) + for (i in seq_along(params)) { + data <- rv$trajectory_stats_df[[params[i]]][!is.na(rv$trajectory_stats_df[[params[i]]])] + if (length(data) >= 3) { + dens <- density(data) + plot(dens, main = paste("Density:", labels[i]), xlab = 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)) + }) + + # Histograms + output$histograms <- renderPlot({ + req(rv$trajectory_stats_df) + + 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") + + par(mfrow = c(2, 3)) + for (i in seq_along(params)) { + data <- rv$trajectory_stats_df[[params[i]]][!is.na(rv$trajectory_stats_df[[params[i]]])] + if (length(data) >= 3) { + hist(data, probability = TRUE, main = paste("Histogram:", labels[i]), + xlab = labels[i], col = "lightgray", border = "darkgray") + lines(density(data), col = "red", lwd = 2) + } + } + par(mfrow = c(1, 1)) + }) + + # Interpretation text + output$interpretation_text <- renderText({ + req(rv$trajectory_stats_df) + + df <- rv$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 + }) +} + +# Run the application +shinyApp(ui = ui, server = server) diff --git a/src/main.Rmd b/src/main.Rmd index 5b68f02..b35c46d 100644 --- a/src/main.Rmd +++ b/src/main.Rmd @@ -43,7 +43,6 @@ 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"]])) -# can get tracks for up to 30 days in the past response <-makeAuthenticatedRequest('tracks/all',query, creds) track_data <- fromJSON(content(response, as = "text", encoding = "UTF-8")) if (!is.null(track_data$path) && length(track_data$path) > 0) { @@ -66,29 +65,383 @@ if (!is.null(track_data$path) && length(track_data$path) > 0) { } ``` -# GUI selection -```{r gui} -icaos <- lapply(departures, function(x) x[["ICAO24"]]) -options <- unlist(icaos) # tcltk needs a character vector - -# Create a GUI list selection -listSelect <- function(options){ - selected_option <- NULL - tryCatch({ - selected_option <- select.list( - title = "Select an aircraft", - choices = options, - preselect = NULL, - multiple = FALSE, - graphics = TRUE - ) - }, error = function(w) { - message('No GUI available') - } +# Trajectory Characteristics Analysis +```{r trajectory-analysis} +if (exists("route_df") && nrow(route_df) > 1) { + + # Convert lat/lon to approximate meters (using simple equirectangular projection) + # Reference point: first coordinate + lat_ref <- route_df$lat[1] + lon_ref <- route_df$lon[1] + + # Convert to meters (approximate) + 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] + + # Create trajr trajectory object + trj <- TrajFromCoords( + data.frame(x = x_meters, y = y_meters, time = time_seconds), + xCol = "x", yCol = "y", timeCol = "time" ) - if (nzchar(selected_option)){ - return(selected_option) - } - return(options[1]) + + # Calculate trajectory characteristics + + # 1. Duration of travel (seconds) + duration <- TrajDuration(trj) + + # 2. Total path length (meters) + path_length <- TrajLength(trj) + + # 3. Diffusion distance (net displacement - straight line from start to end) + diffusion_distance <- TrajDistance(trj) + + # 4. Straightness index (ratio of net displacement to path length, 0-1) + straightness <- TrajStraightness(trj) + + # 5. Mean travel velocity (meters/second) + 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)) + + TrajFractalDimension(trj, stepSizes = step_sizes) + }, error = function(e) { + message("Could not calculate fractal dimension: ", e$message) + NA + }) + + # Create summary data frame + trajectory_characteristics <- data.frame( + Parameter = c( + "Duration of Travel (s)", + "Duration of Travel (min)", + "Path Length (m)", + "Path Length (km)", + "Diffusion Distance (m)", + "Diffusion Distance (km)", + "Straightness Index", + "Mean Travel Velocity (m/s)", + "Mean Travel Velocity (km/h)", + "Fractal Dimension" + ), + Value = c( + round(duration, 2), + round(duration / 60, 2), + round(path_length, 2), + round(path_length / 1000, 2), + round(diffusion_distance, 2), + round(diffusion_distance / 1000, 2), + round(straightness, 4), + round(mean_velocity, 2), + round(mean_velocity * 3.6, 2), + round(fractal_dim, 4) + ) + ) + + print(trajectory_characteristics) + + # Visualize the trajectory using trajr + plot(trj, main = paste("Trajectory of", icao)) + +} else { + message("No valid trajectory data available for analysis") } ``` + +# Statistical Analysis of Multiple Trajectories +```{r multi-trajectory-analysis} +# Function to calculate trajectory characteristics for a single flight +calculate_trajectory_params <- function(icao24, departure_time, creds) { + tryCatch({ + query <- list(icao24 = icao24, time = as.numeric(departure_time)) + response <- makeAuthenticatedRequest('tracks/all', query, creds) + + # Check for HTTP errors + 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) < 2) { + return(NULL) + } + + route_df <- as.data.frame(track_data$path) + colnames(route_df) <- c("time", "lat", "lon", "alt", "heading", "on_ground") + + if (nrow(route_df) < 3) return(NULL) + + # Convert to meters + 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" + ) + + # Calculate parameters + duration <- TrajDuration(trj) + path_length <- TrajLength(trj) + diffusion_dist <- TrajDistance(trj) + straight <- TrajStraightness(trj) + mean_vel <- path_length / duration + + # Fractal dimension + 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) + + return(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) { + message("Error processing ", icao24, ": ", e$message) + 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 +```{r statistical-analysis} +if (exists("trajectory_stats_df") && nrow(trajectory_stats_df) >= 2) { + + # Parameters to analyze + params_to_analyze <- c("diffusion_distance_km", "straightness", "duration_min", + "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)] + if (length(x) < 2) return(NULL) + + data.frame( + Parameter = param_name, + N = length(x), + Mean = round(mean(x), 4), + Variance = round(var(x), 4), + Std_Dev = round(sd(x), 4), + Min = round(min(x), 4), + Q1 = round(quantile(x, 0.25), 4), + Median = round(median(x), 4), + Q3 = round(quantile(x, 0.75), 4), + Max = round(max(x), 4), + IQR = round(IQR(x), 4) + ) + } + + # Calculate statistics for each parameter + stats_list <- list() + for (i in seq_along(params_to_analyze)) { + param <- params_to_analyze[i] + label <- param_labels[i] + stats_list[[i]] <- calc_stats(trajectory_stats_df[[param]], label) + } + + # 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)) + + for (i in seq_along(params_to_analyze)) { + param <- params_to_analyze[i] + label <- param_labels[i] + data <- trajectory_stats_df[[param]][!is.na(trajectory_stats_df[[param]])] + + if (length(data) >= 2) { + boxplot(data, + main = paste("Boxplot:", label), + ylab = label, + col = "lightblue", + border = "darkblue") + + # Add mean point + points(1, mean(data), pch = 18, col = "red", cex = 1.5) + } + } + + par(mfrow = c(1, 1)) + + # Density plots for each parameter + par(mfrow = c(2, 3)) + + for (i in seq_along(params_to_analyze)) { + param <- params_to_analyze[i] + label <- param_labels[i] + data <- trajectory_stats_df[[param]][!is.na(trajectory_stats_df[[param]])] + + if (length(data) >= 3) { + dens <- density(data, na.rm = TRUE) + plot(dens, + main = paste("Density:", label), + xlab = label, + ylab = "Density", + col = "darkblue", + lwd = 2) + 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 = 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)) + + # Histogram with density overlay + par(mfrow = c(2, 3)) + + for (i in seq_along(params_to_analyze)) { + param <- params_to_analyze[i] + label <- param_labels[i] + data <- trajectory_stats_df[[param]][!is.na(trajectory_stats_df[[param]])] + + if (length(data) >= 3) { + hist(data, + probability = TRUE, + main = paste("Histogram:", label), + xlab = label, + col = "lightgray", + border = "darkgray") + + # Overlay density curve + lines(density(data), col = "red", lwd = 2) + } + } + + par(mfrow = c(1, 1)) + +} else { + message("Insufficient trajectory data for statistical analysis (need at least 2 trajectories)") +} +``` + +# Interpretation of Results +```{r interpretation} +if (exists("trajectory_stats_df") && nrow(trajectory_stats_df) >= 2) { + + cat("\n========== INTERPRETATION OF TRAJECTORY PARAMETERS ==========\n\n") + + # Diffusion Distance + dd <- trajectory_stats_df$diffusion_distance_km[!is.na(trajectory_stats_df$diffusion_distance_km)] + if (length(dd) >= 2) { + cat("1. DIFFUSION DISTANCE (Net Displacement):\n") + cat(" - Mean:", round(mean(dd), 2), "km\n") + cat(" - This represents the straight-line distance from origin to destination.\n") + cat(" - High variance (", round(var(dd), 2), ") indicates diverse flight distances.\n\n") + } + + # Straightness + st <- trajectory_stats_df$straightness[!is.na(trajectory_stats_df$straightness)] + if (length(st) >= 2) { + cat("2. STRAIGHTNESS INDEX:\n") + cat(" - Mean:", round(mean(st), 4), "(range 0-1, where 1 = perfectly straight)\n") + cat(" - Values close to 1 indicate efficient, direct flight paths.\n") + cat(" - Lower values suggest deviations due to weather, airspace, or routing.\n\n") + } + + # Duration + dur <- trajectory_stats_df$duration_min[!is.na(trajectory_stats_df$duration_min)] + if (length(dur) >= 2) { + cat("3. DURATION OF TRAVEL:\n") + cat(" - Mean:", round(mean(dur), 2), "minutes\n") + cat(" - Range:", round(min(dur), 2), "-", round(max(dur), 2), "minutes\n") + cat(" - IQR:", round(IQR(dur), 2), "minutes (middle 50% of flights)\n\n") + } + + # Velocity + vel <- trajectory_stats_df$mean_velocity_kmh[!is.na(trajectory_stats_df$mean_velocity_kmh)] + if (length(vel) >= 2) { + cat("4. MEAN TRAVEL VELOCITY:\n") + cat(" - Mean:", round(mean(vel), 2), "km/h\n") + cat(" - Typical commercial aircraft cruise: 800-900 km/h\n") + cat(" - Lower values may include taxi, takeoff, and landing phases.\n\n") + } + + # Fractal Dimension + fd <- trajectory_stats_df$fractal_dimension[!is.na(trajectory_stats_df$fractal_dimension)] + if (length(fd) >= 2) { + cat("5. FRACTAL DIMENSION:\n") + cat(" - Mean:", round(mean(fd), 4), "\n") + cat(" - Value of 1.0 = perfectly straight line\n") + cat(" - Values closer to 2.0 = more complex, space-filling paths\n") + cat(" - Aircraft typically show low fractal dimension (efficient paths).\n\n") + } + + cat("========== END OF ANALYSIS ==========\n") +} +``` +