diff --git a/DESCRIPTION b/DESCRIPTION index 16fc06f..3e1c96b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,6 +7,8 @@ Authors@R: c( person("Robin", "Lovelace", role = c("aut", "ctb"), comment = c(ORCID = "0000-0001-5679-6536")) ) +Author: Zhao Wang [aut, cre], Robin Lovelace [aut, ctb] +Maintainer: Zhao Wang Description: This package provides functions to generate 'core' transport networks from 'complete' transport network datasets. It can take a transport network from 'OpenStreetMap' and saved as an 'sf' object via 'osmextract' or other package @@ -27,6 +29,8 @@ Imports: purrr, units, rlang, + anime, + fmsb Suggests: tibble, tmaptools, diff --git a/R/net_eval.R b/R/net_eval.R index 9894f16..1f7f292 100644 --- a/R/net_eval.R +++ b/R/net_eval.R @@ -23,7 +23,7 @@ compute_spatial_coverage = function(rnet_core, lads, city_name = "City of Edinbu } ### Function: Compute Zone Connectivity -compute_zone_connectivity = function(intermediate_zone, lads, city_name = "City of Edinburgh", buffer_intersection, density_quantile = 0.3) { +compute_zone_connectivity = function(intermediate_zone, lads, city_name = "City of Edinburgh", rnet_core, buffer_distance = 500, density_quantile = 0.3) { # Filter zones by density threshold city_boundary = lads |> filter(LAD23NM == city_name) intermediate_zone = sf::st_intersection(intermediate_zone, city_boundary) @@ -35,8 +35,10 @@ compute_zone_connectivity = function(intermediate_zone, lads, city_name = "City select(InterZone, geometry) |> st_make_valid() - # W = B_union ∩ A_city (already computed outside) - W = buffer_intersection + rnet_core_zone = sf::st_intersection(rnet_core, city_boundary) + network_buffer = st_buffer(rnet_core_zone, dist = buffer_distance) + buffer_union = st_union(network_buffer) + W = sf::st_intersection(buffer_union, city_boundary) # Compute intersections W_i zones = zones |> @@ -84,7 +86,8 @@ compute_zone_connectivity = function(intermediate_zone, lads, city_name = "City graph = g, all_connected = (comp$no == 1), fraction_connected = fraction_connected, - adj_matrix = adj_matrix + adj_matrix = adj_matrix, + buffered_area = W )) } @@ -194,156 +197,343 @@ compute_od_accessibility = function(od_data, rnet_core, lads, city_name = "City )) } -### Function: Compute Directness and Efficiency -compute_directness_efficiency = function(rnet_core, lads, city_name = "City of Edinburgh") { - city_boundary = lads |> filter(LAD23NM == city_name) +### Function: Compute OD-Based Directness and Efficiency +compute_directness_efficiency = function(rnet_core, od_points, lads, city_name = "City of Edinburgh", use_all_od_points = TRUE) { + city_boundary = lads |> filter(LAD23NM == city_name) rnet_core_zone = sf::st_intersection(rnet_core, city_boundary) |> st_cast("LINESTRING") + od_points_zone = od_points[city_boundary, ] + + cat("Network edges:", nrow(rnet_core_zone), "Total OD points:", nrow(od_points_zone), "\n") + + # Check if we have data + if (nrow(rnet_core_zone) == 0 || nrow(od_points_zone) == 0) { + return(list(D = NA, E_glob = NA, E_loc = NA, total_od_count = nrow(od_points_zone), routable_pairs = 0, total_pairs = 0)) + } network_sfnet = as_sfnetwork(rnet_core_zone, directed = FALSE) + + edges_sf = network_sfnet |> activate("edges") |> st_as_sf() + edge_lengths = as.numeric(st_length(edges_sf)) + network_sfnet = network_sfnet |> activate("edges") |> - mutate(weight = as.numeric(st_length(geometry))) + mutate(weight = edge_lengths) g = network_sfnet |> as.igraph() - comp = components(g) - largest_comp_id = which.max(comp$csize) - lcc_nodes = V(g)[comp$membership == largest_comp_id] - g_lcc = induced_subgraph(g, lcc_nodes) - - dist_matrix_g = distances(g_lcc, weights = E(g_lcc)$weight) - nodes_sf = network_sfnet |> activate("nodes") |> st_as_sf() - coords = st_coordinates(nodes_sf) - - node_coords = nodes_sf |> - as.data.frame() |> - bind_cols(data.frame(x = coords[,1], y = coords[,2])) - - node_coords_lcc = node_coords[as.numeric(lcc_nodes), ] - - dist_matrix_e = as.matrix(dist(node_coords_lcc[, c("x","y")])) - - valid_pairs = upper.tri(dist_matrix_g, diag = FALSE) & is.finite(dist_matrix_g) & dist_matrix_g > 0 + # Use ALL OD points for fair comparison across networks + # This ensures same denominator for all network comparisons + od_points_filtered = od_points_zone + n_od = nrow(od_points_filtered) + + cat("Using ALL", n_od, "OD points for fair comparison across networks\n") + + # Find nearest network nodes for all OD points + od_nearest_nodes = st_nearest_feature(od_points_filtered, nodes_sf) + + # Create OD pairs (all combinations of OD points) + od_coords = st_coordinates(od_points_filtered) + + # Calculate city area to scale the number of OD pairs + city_area_km2 = as.numeric(st_area(city_boundary)) / 1e6 # Convert to km² + cat("City area:", round(city_area_km2, 1), "km²\n") + + # Scale max_pairs based on city area + # Base: 10,000 pairs for ~200 km² city (roughly the size of a medium city) + # Scale proportionally with city area + base_area_km2 = 200 + base_max_pairs = 10000 + max_pairs = max(5000, min(50000, round(base_max_pairs * (city_area_km2 / base_area_km2)))) + + cat("Scaled max pairs for this city size:", max_pairs, "\n") + + # Sample OD pairs to make computation manageable + # Use city name as seed to ensure same sampling within same city for fair network comparison + set.seed(abs(sum(utf8ToInt(city_name)))) # Reproducible seed based on city name + + if (n_od > sqrt(max_pairs)) { + sample_size = min(n_od, as.integer(sqrt(max_pairs))) + sampled_indices = sample(1:n_od, sample_size) + od_coords_sample = od_coords[sampled_indices, ] + od_nearest_nodes_sample = od_nearest_nodes[sampled_indices] + od_weights_sample = od_points_filtered$all[sampled_indices] + } else { + od_coords_sample = od_coords + od_nearest_nodes_sample = od_nearest_nodes + od_weights_sample = od_points_filtered$all + } - # Directness - directness_ratios = dist_matrix_e[valid_pairs] / dist_matrix_g[valid_pairs] - D = mean(directness_ratios, na.rm = TRUE) + n_sample = nrow(od_coords_sample) + + # Create all pairs from sampled OD points + od_pairs = expand.grid(origin = 1:n_sample, destination = 1:n_sample) + od_pairs = od_pairs[od_pairs$origin != od_pairs$destination, ] # Remove same-point pairs + + total_pairs = nrow(od_pairs) + cat("Processing", total_pairs, "OD pairs from", n_sample, "sampled OD points\n") + cat("City area-based scaling: area =", round(city_area_km2, 1), "km², max_pairs =", max_pairs, "\n") + cat("Using city-based seed for reproducible sampling within", city_name, "\n") + + # Calculate euclidean distances for OD pairs + origin_coords = od_coords_sample[od_pairs$origin, ] + dest_coords = od_coords_sample[od_pairs$destination, ] + od_euclidean_dist = sqrt(rowSums((origin_coords - dest_coords)^2)) + + # Calculate network distances for OD pairs + origin_nodes = od_nearest_nodes_sample[od_pairs$origin] + dest_nodes = od_nearest_nodes_sample[od_pairs$destination] + + od_network_dist = numeric(length(origin_nodes)) + + for (i in seq_along(origin_nodes)) { + if (origin_nodes[i] != dest_nodes[i]) { + tryCatch({ + path_dist = distances(g, v = origin_nodes[i], to = dest_nodes[i], weights = E(g)$weight) + od_network_dist[i] = path_dist[1, 1] + }, error = function(e) { + od_network_dist[i] = Inf + }) + } else { + od_network_dist[i] = 0 + } + } - # Global Efficiency - sum_inv_dG = sum(1/dist_matrix_g[valid_pairs], na.rm = TRUE) - sum_inv_dE = sum(1/dist_matrix_e[valid_pairs], na.rm = TRUE) - E_glob = sum_inv_dG / sum_inv_dE + # Count successful routing attempts + routable_pairs = sum(is.finite(od_network_dist) & od_euclidean_dist > 0 & od_network_dist > 0) + cat("Successfully routable pairs:", routable_pairs, "out of", total_pairs, + "(", round(100 * routable_pairs / total_pairs, 1), "%)\n") + + # Calculate metrics for ALL pairs (assign 0 efficiency to unsuccessful ones) + # Get weights for all pairs (average of origin and destination weights) + origin_weights = od_weights_sample[od_pairs$origin] + dest_weights = od_weights_sample[od_pairs$destination] + pair_weights = (origin_weights + dest_weights) / 2 + + # Initialize efficiency vectors with zeros + directness_values = numeric(total_pairs) + global_efficiency_values = numeric(total_pairs) + + # Calculate values only for valid pairs + valid_pairs = is.finite(od_network_dist) & od_euclidean_dist > 0 & od_network_dist > 0 + + if (sum(valid_pairs) > 0) { + # OD-based Directness + directness_values[valid_pairs] = od_euclidean_dist[valid_pairs] / od_network_dist[valid_pairs] + + # OD-based Global Efficiency + global_efficiency_values[valid_pairs] = (od_network_dist[valid_pairs]) / (od_euclidean_dist[valid_pairs]) + } - # Local Efficiency - calc_local_eff = function(node_id, distG, distE, g_graph) { - neighbors_idx = as.numeric(neighbors(g_graph, node_id)) - if (length(neighbors_idx) < 2) { - return(NA) + # Calculate weighted averages (unsuccessful pairs contribute 0) + D = weighted.mean(directness_values, pair_weights, na.rm = TRUE) + E_glob = weighted.mean(global_efficiency_values, pair_weights, na.rm = TRUE) + + # OD-based Local Efficiency + # For each OD point, calculate efficiency within its local neighborhood + calc_local_eff_od = function(od_idx, od_coords_sample, od_nearest_nodes_sample, od_weights_sample, g, radius = 2000) { + center_coord = od_coords_sample[od_idx, ] + + # Find OD points within radius + distances_to_center = sqrt(rowSums((sweep(od_coords_sample, 2, center_coord))^2)) + local_indices = which(distances_to_center <= radius & distances_to_center > 0) + + if (length(local_indices) < 2) { + return(0) # Return 0 instead of NA for fair comparison + } + + # Calculate efficiency within local area + local_coords = od_coords_sample[local_indices, ] + local_nodes = od_nearest_nodes_sample[local_indices] + local_weights = od_weights_sample[local_indices] + + # Create pairs within local area + n_local = length(local_indices) + local_pairs = expand.grid(1:n_local, 1:n_local) + local_pairs = local_pairs[local_pairs$Var1 != local_pairs$Var2, ] + + if (nrow(local_pairs) == 0) { + return(0) # Return 0 instead of NA + } + + # Calculate distances for local pairs + local_euc = sqrt(rowSums((local_coords[local_pairs$Var1, ] - local_coords[local_pairs$Var2, ])^2)) + + local_net = numeric(nrow(local_pairs)) + for (i in 1:nrow(local_pairs)) { + node1 = local_nodes[local_pairs$Var1[i]] + node2 = local_nodes[local_pairs$Var2[i]] + if (node1 != node2) { + tryCatch({ + path_dist = distances(g, v = node1, to = node2, weights = E(g)$weight) + local_net[i] = path_dist[1, 1] + }, error = function(e) { + local_net[i] = Inf + }) + } else { + local_net[i] = 0 + } } - distG_sub = distG[neighbors_idx, neighbors_idx] - distE_sub = distE[neighbors_idx, neighbors_idx] - valid_pairs_local = upper.tri(distG_sub, diag=FALSE) & is.finite(distG_sub) & distG_sub > 0 - if (sum(valid_pairs_local) == 0) { - return(NA) + + # Calculate efficiency for ALL local pairs (assign 0 to unsuccessful ones) + local_efficiency_values = numeric(nrow(local_pairs)) + valid_local = is.finite(local_net) & local_euc > 0 & local_net > 0 + + if (sum(valid_local) > 0) { + local_efficiency_values[valid_local] = (1/local_net[valid_local]) / (1/local_euc[valid_local]) } - sum_inv_dG_sub = sum(1/distG_sub[valid_pairs_local], na.rm = TRUE) - sum_inv_dE_sub = sum(1/distE_sub[valid_pairs_local], na.rm = TRUE) - E_glob_i = sum_inv_dG_sub / sum_inv_dE_sub - return(E_glob_i) + + # Calculate weighted average (unsuccessful pairs contribute 0) + local_weights_pairs = (local_weights[local_pairs$Var1] + local_weights[local_pairs$Var2]) / 2 + + return(weighted.mean(local_efficiency_values, local_weights_pairs, na.rm = TRUE)) } - E_glob_i_values = sapply(seq_len(vcount(g_lcc)), function(i) calc_local_eff(i, dist_matrix_g, dist_matrix_e, g_lcc)) - E_loc = mean(E_glob_i_values, na.rm = TRUE) + # Calculate local efficiency for each OD point + E_loc_values = sapply(1:n_sample, function(i) { + calc_local_eff_od(i, od_coords_sample, od_nearest_nodes_sample, od_weights_sample, g) + }) + + E_loc = mean(E_loc_values, na.rm = TRUE) return(list( D = D, E_glob = E_glob, - E_loc = E_loc + E_loc = E_loc, + total_od_count = nrow(od_points_zone), + routable_pairs = routable_pairs, + total_pairs = total_pairs, + routing_success_rate = routable_pairs / total_pairs, + city_area_km2 = city_area_km2, + max_pairs_used = max_pairs )) } generate_radar_chart = function(city_name, - rnet_core, lads, intermediate_zone, - rnet_npt, crs_target, od_data, + rnet_core, rnet_existing, lads, intermediate_zone, + rnet_npt, crs_target, od_data, od_points, dist_threshold = 500, buffer_distance = 500, save_path = NULL) { library(fmsb) - # 1. Spatial Coverage - sp_cov_result = compute_spatial_coverage(rnet_core, lads, city_name = city_name, buffer_distance = buffer_distance) - spatial_coverage = sp_cov_result$coverage * 100 # Convert to percentage - print(paste("Spatial coverage: ", spatial_coverage, "%")) - # 2. Zone Connectivity - zone_conn_result = compute_zone_connectivity(intermediate_zone, lads, city_name = city_name, sp_cov_result$buffered_area) - zone_connectivity = zone_conn_result$fraction_connected * 100 - print(paste("Zone connectivity: ", zone_connectivity, "%")) - - # 3. Cycling Potential Coverage - cp_cov_result = compute_cycling_potential_coverage(rnet_npt, lads, city_name = city_name, rnet_core, crs_target) - cycling_potential_coverage = cp_cov_result$coverage * 100 - citywide_density_ratio = cp_cov_result$coverage_ratio # Might need different scaling if > 1 - print(paste("Cycling potential coverage: ", cycling_potential_coverage, "%")) - - # 4. Population Coverage - pop_cov = compute_population_coverage(intermediate_zone, lads, city_name = city_name, rnet_core, dist_threshold = dist_threshold) - population_coverage = pop_cov * 100 - print(paste("Population coverage: ", population_coverage, "%")) - - # 5. O-D Accessibility - od_result = compute_od_accessibility(od_data, rnet_core, lads, city_name = city_name, dist_threshold = dist_threshold) - od_coverage_count = od_result$od_coverage_count * 100 - od_coverage_demand = od_result$od_coverage_demand * 100 - print(paste("OD coverage (count-based): ", od_coverage_count, "%")) - - # 6. Directness and Efficiency - de_result = compute_directness_efficiency(rnet_core, lads, city_name = city_name) - directness = de_result$D - global_efficiency = de_result$E_glob - local_efficiency = de_result$E_loc - print(paste("Directness: ", directness)) - # Combine metrics - df_metrics_numeric <- data.frame( - SpatialCoverage = spatial_coverage, - ZoneConnectivity = zone_connectivity, - CyclingPotentialCoverage = cycling_potential_coverage, - CitywideDensityRatio = citywide_density_ratio, - PopulationCoverage = population_coverage, - AvgODDistance = od_result$avg_distance, - ODCoverageCountBased = od_coverage_count, - ODCoverageDemandWeighted = od_coverage_demand, - Directness = directness, - GlobalEfficiency = global_efficiency, - LocalEfficiency = local_efficiency - ) + cat("Generating dual-network radar chart for", city_name, "\n") + + # Helper function to calculate all metrics for a given network + calculate_network_metrics = function(rnet, network_name) { + cat(" Calculating metrics for", network_name, "network...\n") + + # 1. Spatial Coverage + sp_cov_result = compute_spatial_coverage(rnet, lads, city_name = city_name, buffer_distance = buffer_distance) + spatial_coverage = sp_cov_result$coverage * 100 + + # 2. Zone Connectivity + zone_conn_result = compute_zone_connectivity(intermediate_zone, lads, city_name = city_name, rnet, buffer_distance = buffer_distance) + zone_connectivity = zone_conn_result$fraction_connected * 100 + + # 3. Cycling Potential Coverage + cp_cov_result = compute_cycling_potential_coverage(rnet_npt, lads, city_name = city_name, rnet, crs_target) + cycling_potential_coverage = cp_cov_result$coverage * 100 + citywide_density_ratio = cp_cov_result$coverage_ratio + + # 4. Population Coverage + pop_cov = compute_population_coverage(intermediate_zone, lads, city_name = city_name, rnet, dist_threshold = dist_threshold) + population_coverage = pop_cov * 100 + + # 5. O-D Accessibility + od_result = compute_od_accessibility(od_data, rnet, lads, city_name = city_name, dist_threshold = dist_threshold) + od_coverage_count = od_result$od_coverage_count * 100 + od_coverage_demand = od_result$od_coverage_demand * 100 + + # 6. OD-Based Directness and Efficiency + de_result = compute_directness_efficiency(rnet, od_points, lads, city_name = city_name) + directness = de_result$D + global_efficiency = de_result$E_glob + local_efficiency = de_result$E_loc + routable_pairs_pct = (de_result$routable_pairs / de_result$total_pairs) * 100 + + # Return all metrics + return(data.frame( + SpatialCoverage = spatial_coverage, + ZoneConnectivity = zone_connectivity, + CyclingPotentialCoverage = cycling_potential_coverage, + CitywideDensityRatio = citywide_density_ratio, + PopulationCoverage = population_coverage, + AvgODDistance = od_result$avg_distance, + ODCoverageCountBased = od_coverage_count, + ODCoverageDemandWeighted = od_coverage_demand, + RoutablePairs = routable_pairs_pct, + Directness = directness, + GlobalEfficiency = global_efficiency, + LocalEfficiency = local_efficiency + )) + } - # Scale metrics for radar chart - df_metrics_scaled = data.frame( - SpatialCoverage = df_metrics_numeric$SpatialCoverage / 100, - ZoneConnectivity = df_metrics_numeric$ZoneConnectivity / 100, - CyclingPotentialCoverage = df_metrics_numeric$CyclingPotentialCoverage / 100, - CitywideDensityRatio = df_metrics_numeric$CitywideDensityRatio / 5, - PopulationCoverage = df_metrics_numeric$PopulationCoverage / 100, - AvgODDistance = df_metrics_numeric$AvgODDistance / 1000, - ODCoverageCountBased = df_metrics_numeric$ODCoverageCountBased / 100, - ODCoverageDemandWeighted = df_metrics_numeric$ODCoverageDemandWeighted / 100, - Directness = df_metrics_numeric$Directness, - GlobalEfficiency = df_metrics_numeric$GlobalEfficiency, - LocalEfficiency = df_metrics_numeric$LocalEfficiency - ) + # Calculate metrics for both networks + core_metrics = calculate_network_metrics(rnet_core, "core") + existing_metrics = calculate_network_metrics(rnet_existing, "existing") + + # Calculate relative performance ratios + directness_ratio = ifelse(existing_metrics$Directness > 0, + core_metrics$Directness / existing_metrics$Directness, + NA) + global_eff_ratio = ifelse(existing_metrics$GlobalEfficiency > 0, + core_metrics$GlobalEfficiency / existing_metrics$GlobalEfficiency, + NA) + local_eff_ratio = ifelse(existing_metrics$LocalEfficiency > 0, + core_metrics$LocalEfficiency / existing_metrics$LocalEfficiency, + NA) + routable_pairs_ratio = ifelse(existing_metrics$RoutablePairs > 0, + core_metrics$RoutablePairs / existing_metrics$RoutablePairs, + NA) + + # Print comparison + cat(" Core network metrics:\n") + cat(" Spatial Coverage:", round(core_metrics$SpatialCoverage, 1), "%\n") + cat(" Routable Pairs:", round(core_metrics$RoutablePairs, 1), "%\n") + cat(" Directness:", round(core_metrics$Directness, 4), "\n") + cat(" Global Efficiency:", round(core_metrics$GlobalEfficiency, 4), "\n") + + cat(" Existing network metrics:\n") + cat(" Spatial Coverage:", round(existing_metrics$SpatialCoverage, 1), "%\n") + cat(" Routable Pairs:", round(existing_metrics$RoutablePairs, 1), "%\n") + cat(" Directness:", round(existing_metrics$Directness, 4), "\n") + cat(" Global Efficiency:", round(existing_metrics$GlobalEfficiency, 4), "\n") + + # Scale metrics for radar chart (excluding efficiency metrics and routable pairs) + scale_metrics = function(df_metrics) { + # Scale and invert Average O-D Distance so that lower distances (better) result in higher values + scaled_distance = df_metrics$AvgODDistance / 1000 # Convert to km + inverted_distance = 1 - pmin(scaled_distance, 1) # Invert and cap at 1 to avoid negative values + + data.frame( + SpatialCoverage = df_metrics$SpatialCoverage / 100, + ZoneConnectivity = df_metrics$ZoneConnectivity / 100, + CyclingPotentialCoverage = df_metrics$CyclingPotentialCoverage / 100, + CitywideDensityRatio = df_metrics$CitywideDensityRatio / 5, + PopulationCoverage = df_metrics$PopulationCoverage / 100, + AvgODDistance = inverted_distance, + ODCoverageCountBased = df_metrics$ODCoverageCountBased / 100, + ODCoverageDemandWeighted = df_metrics$ODCoverageDemandWeighted / 100 + ) + } + + core_scaled = scale_metrics(core_metrics) + cycle_scaled = scale_metrics(existing_metrics) - # Prepare radar chart data - max_vals = rep(1, ncol(df_metrics_scaled)) - min_vals = rep(0, ncol(df_metrics_scaled)) + # Prepare radar chart data with max/min bounds + max_vals = rep(1, ncol(core_scaled)) + min_vals = rep(0, ncol(core_scaled)) + # Combine data for radar chart (max, min, core, cycle) df_radar = rbind( max_vals, min_vals, - df_metrics_scaled + core_scaled, + cycle_scaled ) # Rename columns for radar chart @@ -353,41 +543,112 @@ generate_radar_chart = function(city_name, "Cycling Potential\nCoverage", "Density\nRatio", "Pop.\nCoverage", - "Avg.\nOD\nDistance", + "OD\nAccessibility", "OD\nCoverage\n(Count)", - "OD\nCoverage\n(Demand)", - "Directness", - "Global\nEfficiency", - "Local\nEfficiency" + "OD\nCoverage\n(Demand)" ) # Generate radar chart if (!is.null(save_path)) { - # Open a PNG device to save the plot - png(filename = save_path, width = 800, height = 800) + png(filename = save_path, width = 1200, height = 900) } + # Create the radar chart with both networks radarchart( df_radar, axistype = 1, seg = 5, - pcol = rgb(0.2, 0.5, 0.5, 0.9), - pfcol = rgb(0.2, 0.5, 0.5, 0.5), - plwd = 2, + # Core network in red, cycle network in blue + pcol = c("red", "blue"), + pfcol = c(rgb(1, 0, 0, 0.3), rgb(0, 0, 1, 0.3)), # Semi-transparent fills + plwd = 3, cglcol = "grey", cglty = 1, axislabcol = "grey", - caxislabels = seq(0, max(max_vals), length.out = 6), - title = paste(city_name), - # Increase font sizes - cex.axis = 2, # Axis label size - cex.lab = 2, # Axis title size (if applicable) - cex.main = 2, # Main title size - vlcex = 1.5 + caxislabels = seq(0, 1, length.out = 6), + title = paste(city_name, "- Core (Red) vs Cycle (Blue) Networks"), + cex.axis = 1.2, + cex.main = 1.5, + vlcex = 1.2 ) + # Add legend + legend( + x = "topright", + legend = c("Core Network", "Cycle Network"), + col = c("red", "blue"), + lty = 1:2, + bty = "n" + ) + + # Create performance comparison dataframe + performance_data = data.frame( + Metric = character(), + Core_Value = numeric(), + Existing_Value = numeric(), + Ratio = numeric(), + Performance_Text = character(), + stringsAsFactors = FALSE + ) + + if (!is.na(directness_ratio)) { + performance_data = rbind(performance_data, data.frame( + Metric = "Directness", + Core_Value = core_metrics$Directness, + Existing_Value = existing_metrics$Directness, + Ratio = directness_ratio, + Performance_Text = paste("Core network has", round(directness_ratio, 2), "x better directness") + )) + } + + if (!is.na(global_eff_ratio)) { + performance_data = rbind(performance_data, data.frame( + Metric = "Global_Efficiency", + Core_Value = core_metrics$GlobalEfficiency, + Existing_Value = existing_metrics$GlobalEfficiency, + Ratio = global_eff_ratio, + Performance_Text = paste("Core network has", round(global_eff_ratio, 2), "x better global efficiency") + )) + } + + if (!is.na(local_eff_ratio)) { + performance_data = rbind(performance_data, data.frame( + Metric = "Local_Efficiency", + Core_Value = core_metrics$LocalEfficiency, + Existing_Value = existing_metrics$LocalEfficiency, + Ratio = local_eff_ratio, + Performance_Text = paste("Core network has", round(local_eff_ratio, 2), "x better local efficiency") + )) + } + + if (!is.na(routable_pairs_ratio)) { + performance_data = rbind(performance_data, data.frame( + Metric = "Routable_Pairs", + Core_Value = core_metrics$RoutablePairs, + Existing_Value = existing_metrics$RoutablePairs, + Ratio = routable_pairs_ratio, + Performance_Text = paste("Core network has", round(routable_pairs_ratio, 2), "x better routable pairs") + )) + } + if (!is.null(save_path)) { - dev.off() # Close the PNG device + dev.off() + cat(" Saved radar chart to:", save_path, "\n") } + + # Print performance comparison text + cat(" Performance Comparisons:\n") + for (i in 1:nrow(performance_data)) { + cat(" ", performance_data$Performance_Text[i], "\n") + } + + # Return both sets of metrics and performance comparison dataframe + return(list( + core_metrics = core_metrics, + existing_metrics = existing_metrics, + core_scaled = core_scaled, + cycle_scaled = cycle_scaled, + performance_comparison = performance_data + )) } diff --git a/man/calculate_paths_from_point_dist.Rd b/man/calculate_paths_from_point_dist.Rd index 34d228f..00d70e7 100644 --- a/man/calculate_paths_from_point_dist.Rd +++ b/man/calculate_paths_from_point_dist.Rd @@ -11,7 +11,8 @@ calculate_paths_from_point_dist( maxDistPts = 1500, centroids, path_type = "shortest", - max_path_weight = 10 + max_path_weight = 10, + crs_transform = 27700 ) } \arguments{ diff --git a/network_eval.qmd b/network_eval.qmd index 7090581..b7ad4a7 100644 --- a/network_eval.qmd +++ b/network_eval.qmd @@ -15,15 +15,77 @@ library(tidygraph) library(sfnetworks) library(units) source("R/net_eval.R") +library(osmactive) ``` ```{r} -rnet_core = sf::st_read("network_eval_file/combined_CN_4_2024-12-01_OS.geojson") +check_and_download_file <- function(file_path, download_url, description = "file") { + if (!file.exists(file_path)) { + cat("File", file_path, "does not exist. Downloading", description, "...\n") + + dir.create(dirname(file_path), recursive = TRUE, showWarnings = FALSE) + + tryCatch({ + download.file(download_url, file_path, mode = "wb") + cat("Successfully downloaded", description, "to", file_path, "\n") + }, error = function(e) { + stop("Failed to download ", description, ": ", e$message) + }) + } else { + cat("File", file_path, "already exists. Using existing file.\n") + } +} + +base_url <- "https://github.com/nptscot/corenet/releases/download/net_eval/" + +files_to_check <- list( + list( + path = "network_eval_file/combined_CN_4_2025-05-01_OS.geojson", + url = paste0(base_url, "combined_CN_4_2025-05-01_OS.geojson"), + description = "Core Network GeoJSON" + ), + list( + path = "network_eval_file/SG_IntermediateZone_Bdry_2011.gpkg", + url = paste0(base_url, "SG_IntermediateZone_Bdry_2011.gpkg"), + description = "Intermediate Zone Boundaries" + ), + list( + path = "network_eval_file/combined_od_commute_subset_2025-05-01.Rds", + url = paste0(base_url, "od_izo_sf.Rds"), + description = "OD Data" + ), + list( + path = "network_eval_file/la_regions_scotland_bfe_simplified_2023.geojson", + url = paste0(base_url, "la_regions_scotland_bfe_simplified_2023.geojson"), + description = "Local Authority Boundaries" + ), + list( + path = "network_eval_file/combined_network.gpkg", + url = paste0(base_url, "combined_network.gpkg"), + description = "NPT Combined Network" + ), + list( + path = "network_eval_file/existing_cycle_network.gpkg", + url = paste0(base_url, "existing_cycle_network.gpkg"), + description = "Existing Cycle Network" + ) +) + +# Check and download all required files +for (file_info in files_to_check) { + check_and_download_file(file_info$path, file_info$url, file_info$description) +} +``` + +```{r} +# Now read the files (they should exist at this point) +rnet_core = sf::st_read("network_eval_file/combined_CN_4_2025-05-01_OS.geojson") intermediate_zone = sf::st_read("network_eval_file/SG_IntermediateZone_Bdry_2011.gpkg") intermediate_zone$geometry = intermediate_zone$geom -mapview(intermediate_zone) -od_data = readRDS("network_eval_file/od_izo_sf.Rds") +# mapview(intermediate_zone) +od_data = readRDS("network_eval_file/combined_od_commute_subset_2025-05-01.Rds") +dim(od_data) lads = sf::read_sf("network_eval_file/la_regions_scotland_bfe_simplified_2023.geojson") |> st_transform(27700) @@ -38,24 +100,336 @@ rnet_core = st_transform(rnet_core, crs = crs_target) intermediate_zone = st_transform(intermediate_zone, crs = crs_target) od_data = st_transform(od_data, crs = crs_target) lads = st_transform(lads, crs = crs_target) + +od_points = od_data |> + st_cast("POINT") + +rnet_existing = sf::st_read("network_eval_file/existing_cycle_network.gpkg") +``` + +## Compare original vs jittered OD data for Edinburgh + +```{r} +cat("Loading ready jittered OD dataset...\n") + +od_data_original = readRDS("network_eval_file/od_izo_sf.Rds") |> st_cast("POINT") |> st_transform(27700) +od_jittered = readRDS("network_eval_file/od_commute_subset.Rds") |> st_cast("POINT") |> st_transform(27700) + +edinburgh_city = "City of Edinburgh" + +city_boundary = lads |> filter(LAD23NM == edinburgh_city) + +od_data_original_city = od_data_original[city_boundary, ] +od_jittered_city = od_jittered[city_boundary, ] + +cat("Original OD data dimensions:", nrow(od_data_original_city), "rows,", ncol(od_data_original_city), "columns\n") +cat("Jittered OD data dimensions:", nrow(od_jittered_city), "rows,", ncol(od_jittered_city), "columns\n") + +mapview(od_data_original_city, color = "red") + mapview(od_jittered_city, color = "blue") ``` + + ```{r} city_names = c("Glasgow City", "City of Edinburgh", "Aberdeen City", "Dundee City") +# Generate dual-network radar charts (Core=Red, Cycle=Blue) for (city_name in city_names) { generate_radar_chart( city_name = city_name, rnet_core = rnet_core, + rnet_existing = rnet_existing, lads = lads, intermediate_zone = intermediate_zone, rnet_npt = rnet_npt, crs_target = crs_target, od_data = od_data, + od_points = od_points, dist_threshold = 500, buffer_distance = 500, - save_path = glue::glue("./plot/", city_name , "_radar_chart.jpg") + save_path = glue::glue("./plot/", city_name , "_dual_network_radar_chart.jpg") ) } +``` +Code to generate the existing cycle network (existing_cycle_network.gpkg) +```{r} +# city_names = c("Glasgow City", "City of Edinburgh", "Aberdeen City", "Dundee City") + +# # Initialize an empty list to store cycle networks for each city +# cycle_nets_list = list() + +# # Loop through each city to create cycle networks +# for (i in seq_along(city_names)) { +# city_name = city_names[i] + +# cat("Processing", city_name, "...\n") + +# # Get travel network for the city +# osm = osmactive::get_travel_network("Scotland", boundary = lads[lads$LAD23NM == city_name, ], boundary_type = "clipsrc") +# cycle_net = osmactive::get_cycling_network(osm) +# drive_net = osmactive::get_driving_network(osm) +# cycle_net = osmactive::distance_to_road(cycle_net, drive_net) +# cycle_net = osmactive::classify_cycle_infrastructure(cycle_net) |> st_transform(27700) + +# # Add city identifier to the cycle network +# cycle_net$city = city_name + +# # Store in list +# cycle_nets_list[[i]] = cycle_net + +# cat("Completed", city_name, "\n") +# } + +# # Combine all cycle networks into a single sf object +# rnet_cycle = do.call(rbind, cycle_nets_list) +# # Display summary +# cat("Combined cycle network created with", nrow(rnet_cycle), "features across", length(city_names), "cities\n") + +# # Optional: View the combined network +# mapview::mapview(rnet_cycle, zcol = "city") + +# sf::st_write(rnet_cycle, "network_eval_file/existing_cycle_network.gpkg") ``` + +## OD-Based Efficiency Analysis + +```{r} +source("R/net_eval.R") + +# Initialize results list +results_list = list() + +# Run analysis for all cities +for (i in seq_along(city_names)) { + test_city = city_names[i] + cat("Testing city:", test_city, "\n") + + city_boundary = lads |> filter(LAD23NM == test_city) + od_points_city = od_points[city_boundary, ] + rnet_core_city = sf::st_intersection(rnet_core, city_boundary) + rnet_cycle_city = rnet_existing[city_boundary, ] + + cat("=== TRULY Fair Network Comparison for", test_city, "===\n") + cat("Key principle: Same OD points for all networks, unsuccessful routing = 0 efficiency\n\n") + + # Test with core network + cat("--- Core Network ---\n") + core_results = compute_directness_efficiency(rnet_core_city, od_points_city, lads, city_name = test_city) + + # Test with cycle network + cat("\n--- Cycle Network ---\n") + cycle_results = compute_directness_efficiency(rnet_cycle_city, od_points_city, lads, city_name = test_city) + + # Store results for this city + results_list[[test_city]] = list( + core_network = list( + E_glob = core_results$E_glob, + E_loc = core_results$E_loc, + D = core_results$D, + routing_success_rate = core_results$routing_success_rate, + routable_pairs = core_results$routable_pairs, + total_pairs = core_results$total_pairs + ), + cycle_network = list( + E_glob = cycle_results$E_glob, + E_loc = cycle_results$E_loc, + D = cycle_results$D, + routing_success_rate = cycle_results$routing_success_rate, + routable_pairs = cycle_results$routable_pairs, + total_pairs = cycle_results$total_pairs + ) + ) + + # Compare results + cat("\n=== Comparison Results for", test_city, "===\n") + cat("SAME OD points used for both networks:", core_results$total_od_count, "\n") + cat("SAME total pairs evaluated:", core_results$total_pairs, "\n\n") + + cat("Core Network:\n") + cat("- Successfully routable pairs:", core_results$routable_pairs, "/", core_results$total_pairs, + "(", round(100 * core_results$routing_success_rate, 1), "%)\n") + cat("- Directness (including 0s):", round(core_results$D, 4), "\n") + cat("- Global Efficiency (including 0s):", round(core_results$E_glob, 4), "\n") + cat("- Local Efficiency (including 0s):", round(core_results$E_loc, 4), "\n") + + cat("\nCycle Network:\n") + cat("- Successfully routable pairs:", cycle_results$routable_pairs, "/", cycle_results$total_pairs, + "(", round(100 * cycle_results$routing_success_rate, 1), "%)\n") + cat("- Directness (including 0s):", round(cycle_results$D, 4), "\n") + cat("- Global Efficiency (including 0s):", round(cycle_results$E_glob, 4), "\n") + cat("- Local Efficiency (including 0s):", round(cycle_results$E_loc, 4), "\n") + + # Calculate relative performance + cat("\n=== Relative Performance for", test_city, "===\n") + if (core_results$D > cycle_results$D) { + cat("Core network has", round(core_results$D / cycle_results$D, 2), "x better directness\n") + } else { + cat("Cycle network has", round(cycle_results$D / core_results$D, 2), "x better directness\n") + } + + if (core_results$E_glob > cycle_results$E_glob) { + cat("Core network has", round(core_results$E_glob / cycle_results$E_glob, 2), "x better global efficiency\n") + } else { + cat("Cycle network has", round(cycle_results$E_glob / core_results$E_glob, 2), "x better global efficiency\n") + } + + if (core_results$E_loc > cycle_results$E_loc) { + cat("Core network has", round(core_results$E_loc / cycle_results$E_loc, 2), "x better local efficiency\n") + } else { + cat("Cycle network has", round(cycle_results$E_loc / core_results$E_loc, 2), "x better local efficiency\n") + } + + cat("\n" , rep("=", 50), "\n\n") +} + +# Save results to file +saveRDS(results_list, "network_eval_file/all_cities_efficiency_results.Rds") + +# Print summary of all results +cat("=== SUMMARY OF ALL CITIES ===\n") +for (city in names(results_list)) { + cat("\n", city, ":\n") + cat(" Core Network - E_glob:", round(results_list[[city]]$core_network$E_glob, 4), + ", E_loc:", round(results_list[[city]]$core_network$E_loc, 4), + ", D:", round(results_list[[city]]$core_network$D, 4), "\n") + cat(" Cycle Network - E_glob:", round(results_list[[city]]$cycle_network$E_glob, 4), + ", E_loc:", round(results_list[[city]]$cycle_network$E_loc, 4), + ", D:", round(results_list[[city]]$cycle_network$D, 4), "\n") +} +``` + +| City | Metric | Core Network | Existing Network | Improvement Factor | +|------------------|--------------------|--------------|------------------|---------------------------| +| **Glasgow City** | Directness | 0.0296 | 0.0056 | 5.29 × | +| | Global Efficiency | 0.0634 | 0.0014 | 45.29 × | +| | Local Efficiency | 0.2637 | 0.0380 | 6.94 × | +| | Routable Pairs | 3.2% | 0.2% | 16.00 × | +| **Edinburgh** | Directness | 0.1259 | 0.0049 | 25.71 × | +| | Global Efficiency | 0.4068 | 0.0040 | 101.70 × | +| | Local Efficiency | 0.2374 | 0.0418 | 5.68 × | +| | Routable Pairs | 18.2% | 0.3% | 60.67 × | +| **Aberdeen City**| Directness | 0.4894 | 0.0218 | 22.45 × | +| | Global Efficiency | 0.7914 | 0.0042 | 188.43 × | +| | Local Efficiency | 0.6284 | 0.2356 | 2.67 × | +| | Routable Pairs | 59.6% | 0.4% | 149.00 × | +| **Dundee City** | Directness | 0.0937 | 0.0217 | 4.32 × | +| | Global Efficiency | 0.0583 | 0.0043 | 13.56 × | +| | Local Efficiency | 0.3575 | 0.0695 | 5.14 × | +| | Routable Pairs | 3.6% | 0.6% | 6.00 × | + + +```{r} +library(ggplot2) +library(patchwork) # For combining plots side by side + +# Define performance metrics data +performance_metrics <- list( + "Glasgow City" = list( + directness = c(core = 0.0296, existing = 0.0056, factor = 5.29), + global_eff = c(core = 0.0634, existing = 0.0014, factor = 45.29), + local_eff = c(core = 0.2637, existing = 0.0380, factor = 6.94), + routable = c(core = 3.2, existing = 0.2, factor = 16.00) + ), + "City of Edinburgh" = list( + directness = c(core = 0.1259, existing = 0.0049, factor = 25.71), + global_eff = c(core = 0.4068, existing = 0.0040, factor = 101.70), + local_eff = c(core = 0.2374, existing = 0.0418, factor = 5.68), + routable = c(core = 18.2, existing = 0.3, factor = 60.67) + ), + "Aberdeen City" = list( + directness = c(core = 0.4894, existing = 0.0218, factor = 22.45), + global_eff = c(core = 0.7914, existing = 0.0042, factor = 188.43), + local_eff = c(core = 0.6284, existing = 0.2356, factor = 2.67), + routable = c(core = 59.6, existing = 0.4, factor = 149.00) + ), + "Dundee City" = list( + directness = c(core = 0.0937, existing = 0.0217, factor = 4.32), + global_eff = c(core = 0.0583, existing = 0.0043, factor = 13.56), + local_eff = c(core = 0.3575, existing = 0.0695, factor = 5.14), + routable = c(core = 3.6, existing = 0.6, factor = 6.00) + ) +) + +# Create plots for each city +for (city_name in city_names) { + cat("Creating plots for", city_name, "...\n") + + # Filter city boundary + city_boundary = lads |> filter(LAD23NM == city_name) + + # Filter networks for the city + rnet_core_city = sf::st_intersection(rnet_core, city_boundary) + rnet_existing_city = sf::st_intersection(rnet_existing, city_boundary) + + # Get performance metrics for this city + metrics = performance_metrics[[city_name]] + + # Create performance text + perf_text = paste( + sprintf("Directness: %.3f vs %.3f (%.1fx)", metrics$directness["core"], metrics$directness["existing"], metrics$directness["factor"]), + sprintf("Global Eff: %.3f vs %.3f (%.1fx)", metrics$global_eff["core"], metrics$global_eff["existing"], metrics$global_eff["factor"]), + sprintf("Local Eff: %.3f vs %.3f (%.1fx)", metrics$local_eff["core"], metrics$local_eff["existing"], metrics$local_eff["factor"]), + sprintf("Routable: %.1f%% vs %.1f%% (%.1fx)", metrics$routable["core"], metrics$routable["existing"], metrics$routable["factor"]), + sep = "\n" + ) + + # Create core network plot + plot_core <- ggplot() + + geom_sf(data = city_boundary, fill = "lightgray", color = "darkgray", alpha = 0.3) + + geom_sf(data = rnet_core_city, color = "red", size = 0.8) + + ggtitle(paste("Core Network -", city_name)) + + theme_void() + + theme( + plot.title = element_text(hjust = 0.5, size = 12, face = "bold"), + plot.margin = margin(10, 10, 10, 10) + ) + + # Create existing network plot + plot_existing <- ggplot() + + geom_sf(data = city_boundary, fill = "lightgray", color = "darkgray", alpha = 0.3) + + geom_sf(data = rnet_existing_city, color = "blue", size = 0.8) + + ggtitle(paste("Existing Network -", city_name)) + + theme_void() + + theme( + plot.title = element_text(hjust = 0.5, size = 12, face = "bold"), + plot.margin = margin(10, 10, 10, 10) + ) + + # Create performance metrics plot + plot_metrics <- ggplot() + + annotate("text", x = 0.5, y = 0.5, + label = paste("Performance Comparison\n(Core vs Existing)\n\n", perf_text), + hjust = 0.5, vjust = 0.5, size = 3.5, + fontface = "bold", color = "darkblue") + + xlim(0, 1) + ylim(0, 1) + + theme_void() + + theme( + panel.border = element_rect(color = "lightgray", fill = NA, size = 1), + plot.margin = margin(10, 10, 10, 10) + ) + + # Combine plots: networks side by side, metrics below + combined_plot <- (plot_core + plot_existing) / plot_metrics + + plot_layout(heights = c(2, 1)) + + plot_annotation( + title = paste("Network Comparison:", city_name), + theme = theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold")) + ) + + # Display the plot + print(combined_plot) + + # Save the plot + ggsave( + filename = paste0("plot/network_comparison_", gsub(" ", "_", city_name), ".png"), + plot = combined_plot, + width = 12, + height = 9, + dpi = 300, + bg = "white" + ) + + cat("Saved plot for", city_name, "\n\n") +} +``` \ No newline at end of file diff --git a/plot/Aberdeen City_radar_chart.jpg b/plot/Aberdeen City_radar_chart.jpg index 0fdebbb..a9126fa 100644 Binary files a/plot/Aberdeen City_radar_chart.jpg and b/plot/Aberdeen City_radar_chart.jpg differ diff --git a/plot/City of Edinburgh_radar_chart.jpg b/plot/City of Edinburgh_radar_chart.jpg index 67e8154..483d3ab 100644 Binary files a/plot/City of Edinburgh_radar_chart.jpg and b/plot/City of Edinburgh_radar_chart.jpg differ diff --git a/plot/Dundee City_radar_chart.jpg b/plot/Dundee City_radar_chart.jpg index ce3d5e3..2ee0675 100644 Binary files a/plot/Dundee City_radar_chart.jpg and b/plot/Dundee City_radar_chart.jpg differ diff --git a/plot/Glasgow City_radar_chart.jpg b/plot/Glasgow City_radar_chart.jpg index 1dfa001..be70a9f 100644 Binary files a/plot/Glasgow City_radar_chart.jpg and b/plot/Glasgow City_radar_chart.jpg differ diff --git a/plot/Stirling_radar_chart.jpg b/plot/Stirling_radar_chart.jpg deleted file mode 100644 index caf0de7..0000000 Binary files a/plot/Stirling_radar_chart.jpg and /dev/null differ