utility trip checks

.gitignore

@@ -45,3 +45,5 @@ OS_Scotland_Network.geojson
 npt_net_qgis
 scottraffic/
 *.geojson
+
+/.quarto/

R/get_orcp_cn.R

@@ -1,3 +1,9 @@
+library(sfnetworks)
+library(sf)
+library(tidygraph)
+library(tidyverse)
+library(igraph)
+
 orcp_network = function(area, NPT_zones, length_threshold = 10000, percentile_value = 0.6) {
 
   osm = osmactive::get_travel_network("Scotland", boundary = area, boundary_type = "clipsrc")
@@ -99,4 +105,245 @@ orcp_network = function(area, NPT_zones, length_threshold = 10000, percentile_va
     filter(st_intersects(geometry, summarized_data$geometry, sparse = FALSE) |> apply(1, any))
 
     return(cycle_net_NPT_filtered)
-}
+}
+
+find_component= function(rnet, threshold = 50) {
+
+  sf::st_crs(rnet) = 27700
+
+  # Calculate the distance matrix between features
+  dist_matrix = sf::st_distance(rnet)
+
+  # Convert the threshold to units of meters
+  threshold = units::set_units(threshold, "m")
+
+  # Create a connectivity matrix where connections are based on the threshold distance
+  connectivity_matrix = Matrix::Matrix(dist_matrix < threshold, sparse = TRUE)
+
+  # Create an undirected graph from the adjacency matrix
+  graph = igraph::graph_from_adjacency_matrix(connectivity_matrix, mode = "undirected", diag = FALSE)
+
+  # Find the connected components in the graph
+  components = igraph::components(graph)
+
+  # Assign component membership to the road network
+  rnet$component = components$membership
+
+  # Return the updated road network with component membership
+  return(rnet)
+}
+
+find_endpoints = function(rnet) {
+  # Remove duplicate lines in rnet
+  merged_lines = sf::st_union(rnet)
+
+  # Cast to LINESTRING
+  merged_lines = sf::st_cast(merged_lines, "LINESTRING")
+
+  # Create an empty list to store points
+  junctions = list()
+
+  # Extract the first and last coordinates of each line
+  all_points = do.call(rbind,
+    lapply(sf::st_geometry(merged_lines), function(line) {
+      coords = sf::st_coordinates(line)
+      rbind(coords[1, , drop = FALSE], coords[nrow(coords), , drop = FALSE])
+    })
+  )
+
+  if (is.matrix(all_points)) {
+    all_points = as.data.frame(all_points)
+  }
+
+  all_points = all_points[, 1:2]
+  names(all_points) = c("X", "Y")
+
+  # Count occurrences of each point
+  points_counted = all_points |>
+    dplyr::group_by(X, Y) |>
+    dplyr::summarise(Count = dplyr::n(), .groups = 'drop')
+
+  # Filter for unique points (occurring only once)
+  unique_points = points_counted |>
+    dplyr::filter(Count == 1)
+
+  # Convert to sf object
+  points_sf = sf::st_as_sf(unique_points, coords = c("X", "Y"), crs = sf::st_crs(rnet))
+
+  return(points_sf)
+}
+
+
+find_nearest_points = function(points_sf, rnet, segment_length = 10, dist = 100) {
+    rnet = stplanr::line_cast(rnet)
+    rnet = stplanr::line_segment(rnet, segment_length = segment_length)   
+    rnet_points_df = as.data.frame(sf::st_coordinates(rnet))
+    names(rnet_points_df) = c("X", "Y", "L1")
+    rnet_points_sf = st_as_sf(rnet_points_df, coords = c("X", "Y"), crs = st_crs(rnet))
+    points_sf_buffer = st_buffer(points_sf, dist = dist)
+    rnet_points_sf_within_buffer = sf::st_intersection(rnet_points_sf, points_sf_buffer)
+    # Find the nearest point in os_points_sf_within_buffer for each point in points_sf
+    nearest_points = nngeo::st_nn(points_sf, rnet_points_sf_within_buffer, k = 1, returnDist = TRUE)
+
+    # Extract the nearest points and distances
+    nearest_indices = unlist(nearest_points$nn)
+    nearest_distances = unlist(nearest_points$dist)
+
+    # Ensure nearest_indices are within the range of points_sf
+    if (all(nearest_indices <= nrow(rnet_points_sf_within_buffer))) {
+    # Add the index of points_sf to the nearest_os_points
+    nearest_os_points = rnet_points_sf_within_buffer[nearest_indices, ]
+    nearest_os_points$index_points_sf = 1:length(nearest_indices)
+    } else {
+    stop("The nearest_indices are out of bounds. Please check your data.")
+    }
+    return(nearest_os_points)
+}
+
+calculate_paths_from_point_dist = function(network, point, target_point, path_type = "shortest") {
+    # Convert points to sfc
+    point_geom = sf::st_as_sfc(point, crs = st_crs(network))
+    target_geom = sf::st_as_sfc(target_point, crs = st_crs(network))
+
+    # Find nearest network nodes
+    start_node = sf::st_nearest_feature(network, point_geom)
+    end_node = sf::st_nearest_feature(network, target_geom)
+
+    # Ensure network is in the right format and calculate weights
+    network = network |>
+        sf::st_cast("LINESTRING") |>
+        sfnetworks::as_sfnetwork(directed = FALSE) |>
+        sfnetworks::activate("edges") |>
+        dplyr::mutate(weight = sf::st_length(geometry))
+
+    # Calculate paths using weights
+    paths = sfnetworks::st_network_paths(
+        network, 
+        from = start_node, 
+        to = end_node, 
+        weights = "weight",
+        type = path_type
+    )
+
+    # Extract the edge paths
+    edges_in_paths = paths |>
+        dplyr::pull(edge_paths) |>
+        unlist() |>
+        unique()
+
+    # Result as an sf object
+    result = network |>
+        dplyr::slice(edges_in_paths) |>
+        sf::st_as_sf()
+
+    return(result)
+}
+
+compute_shortest_paths = function(points_sf, osm_points, rnet, segment_length = 10) {
+  hull_sf = st_convex_hull(st_combine(osm_points))
+  buffer_distance = 500 
+  buffered_hull_sf = st_buffer(hull_sf, dist = buffer_distance)
+
+  rnet = rnet[sf::st_union(buffered_hull_sf), , op = sf::st_intersects]  
+
+  # Process the OSM data
+  rnet = stplanr::line_cast(rnet)
+  rnet_split = stplanr::line_segment(rnet, segment_length = segment_length)
+
+  # Build the network
+  network = rnet_split |>
+    sf::st_cast("LINESTRING") |>
+    sf::st_transform(27700) |>
+    sfnetworks::as_sfnetwork(directed = FALSE) |>
+    sfnetworks::activate("edges") |>
+    dplyr::mutate(weight = sf::st_length(geometry))
+
+  # Convert sfnetwork to igraph for component analysis
+  network_igraph = as.igraph(network, directed = FALSE)
+  components = igraph::components(network_igraph)
+
+  # Assign component IDs to nodes and edges
+  network = network |>
+    activate("nodes") |>
+    mutate(component_id = as.numeric(components$membership))
+
+  network = network |>
+    activate("edges") |>
+    mutate(component_id = as.numeric(components$membership[head_of(network_igraph, E(network_igraph))]))
+
+  # Initialize an empty list to store paths
+  all_paths = list()
+
+  # Loop through all pairs of points
+  for (i in 1:nrow(points_sf)) {
+    start_point = points_sf[i, ]
+    end_point = osm_points |> filter(index_points_sf == i)
+
+    start_node = sf::st_nearest_feature(start_point, network |> activate("nodes") |> st_as_sf())
+    end_node = sf::st_nearest_feature(end_point, network |> activate("nodes") |> st_as_sf())
+
+    # Compute the shortest path
+    paths_from_point = sfnetworks::st_network_paths(
+      network,
+      from = start_node,
+      to = end_node,
+      weights = "weight",
+      type = "shortest"
+    )
+
+    # Extract edges from the path for visualization
+    edges_in_path = network |>
+      activate("edges") |>
+      filter(row_number() %in% unlist(paths_from_point$edge_paths[[1]])) |>
+      sf::st_as_sf()
+
+    # Store the path information in the list
+    all_paths[[i]] = list(
+      start_point = start_point,
+      end_point = end_point,
+      path_edges = edges_in_path,
+      start_node_geom = network |> activate("nodes") |> slice(start_node) |> st_as_sf(),
+      end_node_geom = network |> activate("nodes") |> slice(end_node) |> st_as_sf()
+    )
+  }
+
+  return(all_paths)
+}
+
+fix_net_connectivity <- function(file_name, input_dir, output_dir, gisBase, gisDbase, location, mapset) {
+  library(rgrass)
+
+  # Initialize GRASS GIS environment
+  initGRASS(gisBase = gisBase, 
+            gisDbase = gisDbase, 
+            location = location, 
+            mapset = mapset, 
+            override = TRUE)
+
+  execGRASS("g.gisenv", flags = "s")
+
+  # Construct input and output file paths
+  input_file <- glue::glue("{input_dir}/{file_name}")
+  output_file <- glue::glue("{output_dir}/{file_name}_fixed.geojson")
+
+  # Import network data from GeoJSON
+  execGRASS("v.in.ogr", flags = c("overwrite", "o"), 
+            parameters = list(input = input_file, 
+                              output = "network"))
+
+  # Clean the network
+  execGRASS("v.clean", flags = c("overwrite"),
+            parameters = list(input = "network", 
+                              output = "fixed_network", 
+                              tool = "break", 
+                              threshold = 10))
+
+  # Export the cleaned network to GeoJSON
+  execGRASS("v.out.ogr", flags = c("overwrite"), 
+            parameters = list(input = "fixed_network", 
+                              output = output_file, 
+                              format = "GeoJSON"))
+
+  cat("Network", output_file, " is cleaned and export completed successfully.\n")
+}
+

_targets.R

@@ -835,11 +835,11 @@ list(
     # geo_code1 and 2 refere to non-Data Zone ids
     end_point = lwgeom::st_endpoint(od_utility_combined)
     end_point = sf::st_join(sf::st_as_sf(end_point), zones)
-    od_utility_combined$endDZ = end_point$DataZone
+    od_utility_combined$geo_code2 = end_point$DataZone
 
     start_point = lwgeom::st_startpoint(od_utility_combined)
     start_point = sf::st_join(sf::st_as_sf(start_point), zones)
-    od_utility_combined$startDZ = start_point$DataZone
+    od_utility_combined$geo_code1 = start_point$DataZone
 
     od_utility_combined
   }),
@@ -975,11 +975,11 @@ list(
   }),
 
   # Utility Zone stats ---------------------------------------------------------
-
   tar_target(utility_stats_baseline, {
-    stats = sf::st_drop_geometry(od_utility_combined)
-     stats = stats[, c(
-      "startDZ", "endDZ", "purpose", "all", "car",
+    stats = sf::st_drop_geometry(uptake_utility_fastest)
+    # Group by start and 
+    stats = stats[, c(
+      "geo_code1", "geo_code2", "purpose", "all", "car",
       "foot", "bicycle", "public_transport", "taxi"
     )]
 
@@ -999,27 +999,27 @@ list(
     stats = rbind(stats_shopping, stats_leisure, stats_visiting)
 
     stats_orig = stats |>
-      dplyr::select(!endDZ) |>
-      dplyr::group_by(startDZ, purpose) |>
+      dplyr::select(!geo_code2) |>
+      dplyr::group_by(geo_code1, purpose) |>
       dplyr::summarise_all(sum, na.rm = TRUE)
 
     stats_dest = stats |>
-      dplyr::select(!startDZ) |>
-      dplyr::group_by(endDZ, purpose) |>
+      dplyr::select(!geo_code1) |>
+      dplyr::group_by(geo_code2, purpose) |>
       dplyr::summarise_all(sum, na.rm = TRUE)
 
     stats_orig$purpose = paste0(stats_orig$purpose, "_orig")
     stats_dest$purpose = paste0(stats_dest$purpose, "_dest")
 
     stats_orig = tidyr::pivot_wider(stats_orig,
-      id_cols = startDZ,
+      id_cols = geo_code1,
       names_from = "purpose",
       values_from = all:taxi,
       names_glue = "{purpose}_{.value}"
     )
 
     stats_dest = tidyr::pivot_wider(stats_dest,
-      id_cols = endDZ,
+      id_cols = geo_code2,
       names_from = "purpose",
       values_from = all:taxi,
       names_glue = "{purpose}_{.value}"
@@ -1030,6 +1030,7 @@ list(
     stats_all = dplyr::full_join(stats_orig, stats_dest, by = "DataZone")
     stats_all
   }),
+
   tar_target(utility_stats_fastest, {
     make_utility_stats(uptake_utility_fastest, "fastest", zones)
   }),