pkgs <- c(“sf”, “ggplot2”, “rnaturalearth”, “rnaturalearthdata”, “igraph”, “dplyr”, “ggspatial”, “osmextract”, “httr”, “jsonlite”)
for (p in pkgs) { if (!requireNamespace(p, quietly = TRUE)) install.packages(p) } library(sf); library(ggplot2) library(rnaturalearth); library(rnaturalearthdata) library(igraph); library(dplyr) library(ggspatial); library(osmextract) library(httr); library(jsonlite)
set.seed(42) utm32 <- st_crs(32632)
denmark_wgs <- ne_countries(country = “Denmark”, scale = “medium”, returnclass = “sf”) |> st_geometry() denmark_utm <- st_transform(denmark_wgs, utm32)
cph_wgs <- st_sfc(st_point(c(12.5683, 55.6761)), crs = 4326) cph_utm <- st_transform(cph_wgs, utm32) |> st_coordinates()
depot <- data.frame(name = “Copenhagen”, X = cph_utm[1,“X”], Y = cph_utm[1,“Y”], depot = TRUE, lon = 12.5683, lat = 55.6761)
message(“Generating 100 towns…”) bbox <- st_bbox(denmark_utm) towns <- data.frame()
while (nrow(towns) < 100) { xs <- runif(3000, bbox[“xmin”], bbox[“xmax”]) ys <- runif(3000, bbox[“ymin”], bbox[“ymax”]) pts <- st_as_sf(data.frame(X = xs, Y = ys), coords = c(“X”,“Y”), crs = utm32) ok <- st_within(pts, st_union(denmark_utm), sparse = FALSE)[,1] cr <- st_coordinates(pts[ok,]) towns <- rbind(towns, data.frame(X = cr[,“X”], Y = cr[,“Y”])) } towns <- towns[1:100,]
towns_sf <- st_as_sf(towns, coords = c(“X”,“Y”), crs = utm32) towns_wgs <- st_transform(towns_sf, 4326) towns_ll <- st_coordinates(towns_wgs)
towns\(name <- paste0("Town_", sprintf("%03d", 1:100)) towns\)depot <- FALSE towns\(lon <- towns_ll[,"X"] towns\)lat <- towns_ll[,“Y”]
all_nodes <- bind_rows(depot, towns) rownames(all_nodes) <- NULL n_nodes <- nrow(all_nodes) # 101 message(sprintf(“Total stops: %d (depot + 100 towns)”, n_nodes))
message(“Computing road-distance matrix via OSRM (chunked)…”) message(” This will take ~2-4 minutes (HTTP calls to router.project-osrm.org)…“)
coords_str <- paste( sprintf(“%.6f,%.6f”, all_nodes\(lon, all_nodes\)lat), collapse = “;” )
osrm_url <- paste0( “http://router.project-osrm.org/table/v1/driving/”, coords_str, “?annotations=distance” )
message(” Sending OSRM table request…“) resp <- tryCatch( GET(osrm_url, timeout(120)), error = function(e) { message(” OSRM request failed: “, conditionMessage(e)) NULL } )
if (!is.null(resp) && status_code(resp) == 200) { osrm_data <- fromJSON(content(resp, “text”, encoding = “UTF-8”)) road_dist_mat <- matrix( unlist(osrm_data$distances), nrow = n_nodes, ncol = n_nodes, byrow = TRUE ) message(” OSRM distance matrix complete.”) } else { # Straight line Euclidean distances in UTM metres message(” OSRM unavailable — using straight-line distances as fallback.”) coords_mat <- as.matrix(all_nodes[, c(“X”,“Y”)]) road_dist_mat <- as.matrix(dist(coords_mat)) }
max_d <- max(road_dist_mat[is.finite(road_dist_mat) & road_dist_mat > 0], na.rm = TRUE) road_dist_mat[!is.finite(road_dist_mat) | is.na(road_dist_mat)] <- max_d * 10 diag(road_dist_mat) <- 0
message(“Solving TSP via MST-DFS…”) g_tsp <- graph_from_adjacency_matrix(road_dist_mat, mode = “undirected”, weighted = TRUE, diag = FALSE) mst_g <- mst(g_tsp) dfs_res <- dfs(mst_g, root = 1, neimode = “all”, unreachable = FALSE) tour <- c(as.integer(dfs_res$order), 1L)
message(“Applying 2-opt improvement…”) two_opt <- function(tour, dm) { n <- length(tour) - 1L improved <- TRUE while (improved) { improved <- FALSE for (i in 1:(n - 2L)) { for (j in (i + 2L):n) { if (i == 1L && j == n) next if ((dm[tour[i], tour[j]] + dm[tour[i+1L], tour[j+1L]]) < (dm[tour[i], tour[i+1L]] + dm[tour[j], tour[j+1L]]) - 1e-6) { tour[(i+1L):j] <- rev(tour[(i+1L):j]) improved <- TRUE } } } } tour } tour <- two_opt(tour, road_dist_mat)
unique_tour <- tour[1:n_nodes] stop_num <- integer(n_nodes) for (i in seq_along(unique_tour)) stop_num[unique_tour[i]] <- i - 1L all_nodes\(stop_num <- stop_num all_nodes\)map_label <- ifelse(all_nodes$depot, “Copenhagen”, as.character(stop_num))
tour_length_km <- round( sum(road_dist_mat[cbind(head(tour,-1), tail(tour,-1))]) / 1000, 1) cat(sprintf(“length: %.1f km”, tour_length_km))
message(“Fetching road geometry for each tour leg from OSRM…”)
get_osrm_route <- function(lon1, lat1, lon2, lat2) { url <- sprintf( “http://router.project-osrm.org/route/v1/driving/%.6f,%.6f;%.6f,%.6f?geometries=geojson&overview=full”, lon1, lat1, lon2, lat2 ) resp <- tryCatch(GET(url, timeout(30)), error = function(e) NULL) if (is.null(resp) || status_code(resp) != 200) return(NULL) parsed <- fromJSON(content(resp, “text”, encoding = “UTF-8”), simplifyVector = FALSE) if (parsed\(code != "Ok") return(NULL) coords <- parsed\)routes[[1]]\(geometry\)coordinates mat <- do.call(rbind, lapply(coords, function(c) c(c[[1]], c[[2]]))) if (nrow(mat) < 2) return(NULL) st_linestring(mat) }
n_legs <- length(tour) - 1 leg_geoms <- vector(“list”, n_legs)
for (i in seq_len(n_legs)) { a <- tour[i]; b <- tour[i + 1] geom <- get_osrm_route( all_nodes\(lon[a], all_nodes\)lat[a], all_nodes\(lon[b], all_nodes\)lat[b] ) if (is.null(geom)) { # straight-line fallback in WGS84 geom <- st_linestring(rbind( c(all_nodes\(lon[a], all_nodes\)lat[a]), c(all_nodes\(lon[b], all_nodes\)lat[b]) )) } leg_geoms[[i]] <- geom if (i %% 20 == 0) message(sprintf(” Leg %d / %d done”, i, n_legs)) }
route_sfc <- st_sfc(leg_geoms, crs = 4326) |> st_transform(utm32) route_sf <- st_sf(geometry = route_sfc) message(” Route geometry complete.”)
message(“Downloading road background layer from Geofabrik…”) options(timeout = 600)
road_display_types <- c(“motorway”,“trunk”,“primary”,“secondary”)
roads_bg <- tryCatch({ raw <- oe_get( place = “Denmark”, layer = “lines”, provider = “geofabrik”, force_download = FALSE, force_vectortranslate = FALSE, query = “SELECT * FROM lines WHERE highway IS NOT NULL”, max_file_size = 5e8, quiet = TRUE ) raw |> filter(highway %in% road_display_types) |> filter(as.character(st_geometry_type(geometry)) == “LINESTRING”) |> filter(!st_is_empty(geometry)) |> st_transform(utm32) |> select(highway, geometry) }, error = function(e) { message(” Road background unavailable: “, conditionMessage(e)) NULL })
towns_df <- all_nodes[!all_nodes\(depot,] depot_df <- all_nodes[ all_nodes\)depot,]
message(“Building map…”)
p <- ggplot() +
# Country – white fill, black border geom_sf(data = denmark_utm, fill = “white”, colour = “black”, linewidth = 0.5) +
# Background roads { if (!is.null(roads_bg)) geom_sf(data = roads_bg, colour = “#bbbbbb”, linewidth = 0.12, alpha = 0.7) else list() } +
# TSP route – red geom_sf(data = route_sf, colour = “#CC0000”, linewidth = 0.55, alpha = 0.92) +
# Town circles – blue geom_point(data = towns_df, aes(x = X, y = Y), shape = 21, fill = “#2171b5”, colour = “black”, size = 2.0, stroke = 0.4) +
# Copenhagen circle – green geom_point(data = depot_df, aes(x = X, y = Y), shape = 21, fill = “#238b45”, colour = “black”, size = 4.5, stroke = 0.9) +
# Town stop numbers geom_text(data = towns_df, aes(x = X, y = Y, label = map_label), size = 1.7, colour = “black”, hjust = 0, nudge_x = 3200, check_overlap = TRUE) +
# Copenhagen label geom_text(data = depot_df, aes(x = X, y = Y, label = map_label), size = 3.2, colour = “black”, fontface = “bold”, nudge_y = 16000, hjust = 0.5) +
annotation_scale(location = “bl”, width_hint = 0.22, unit_category = “metric”) + annotation_north_arrow(location = “tr”, style = north_arrow_fancy_orienteering()) +
labs( title = “COCACOLA OPTIMAL DISTRIBUTION ROUTE IN DENMARK”, subtitle = sprintf( “100 towns | Distribution Centre: Copenhagen | Road distance: %.1f km | UTM Zone 32N (EPSG:32632)”, tour_length_km), x = “Easting”, y = “Northing” ) +
theme_minimal(base_size = 10) + theme( plot.title = element_text(size = 11, face = “bold”, colour = “black”, hjust = 0.5, margin = margin(b = 4), decoration = “underline”), plot.subtitle = element_text(size = 7, colour = “black”, hjust = 0.5), plot.caption = element_blank(), panel.background = element_rect(fill = “white”, colour = NA), plot.background = element_rect(fill = “white”, colour = NA), panel.grid.major = element_line(colour = “#cccccc”, linewidth = 0.25), axis.text = element_text(size = 7), axis.title = element_text(size = 8) )
print(p)
pdf_dir <- “C:/Users/User/Desktop/CLASS NOTES/4.2 Notes/GIS Applications/Assignments/Assignment 5” if (!dir.exists(pdf_dir)) dir.create(pdf_dir, recursive = TRUE) pdf_file <- file.path(pdf_dir, “denmark_tsp_route.pdf”)
pdf(pdf_file, width = 11.69, height = 8.27, paper = “a4r”) print(p) dev.off() cat(sprintf(“saved to:%s”, pdf_file))
cat(“— Tour order —”) for (i in seq_along(unique_tour)) { idx <- unique_tour[i] cat(sprintf(” Stop %3d : %s“, i - 1L, ifelse(all_nodes\(depot[idx], "Copenhagen (Depot)", all_nodes\)name[idx]))) } cat(” –> Return to Copenhagen“)