# Packages
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)
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(rnaturalearth); library(rnaturalearthdata)
## Warning: package 'rnaturalearthdata' was built under R version 4.3.3
##
## Attaching package: 'rnaturalearthdata'
## The following object is masked from 'package:rnaturalearth':
##
## countries110
library(igraph); library(dplyr)
## Warning: package 'igraph' was built under R version 4.3.3
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
## Warning: package 'dplyr' was built under R version 4.3.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:igraph':
##
## as_data_frame, groups, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggspatial); library(osmextract)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
## Check the package website, https://docs.ropensci.org/osmextract/, for more details.
library(httr); library(jsonlite)
## Warning: package 'httr' was built under R version 4.3.3
## Warning: package 'jsonlite' was built under R version 4.3.3
set.seed(42)
utm32 <- st_crs(32632)
# 1.Denmark boundary
denmark_wgs <- ne_countries(country = "Denmark", scale = "medium",
returnclass = "sf") |> st_geometry()
denmark_utm <- st_transform(denmark_wgs, utm32)
# 2. Copenhagen distribution centre
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)
# 100 towns within Denmark
message("Generating 100 towns...")
## 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,]
# Convert UTM towns back to WGS84 lon/lat for OSRM
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"]
# 4. All nodes
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))
## Total stops: 101 (depot + 100 towns)
# 5. Road-distance matrix via OSRM public API
message("Computing road-distance matrix via OSRM (chunked)...")
## Computing road-distance matrix via OSRM (chunked)...
message(" This will take ~2-4 minutes (HTTP calls to router.project-osrm.org)...")
## This will take ~2-4 minutes (HTTP calls to router.project-osrm.org)...
# Build coordinate string: "lon,lat;lon,lat;..."
coords_str <- paste(
sprintf("%.6f,%.6f", all_nodes$lon, all_nodes$lat),
collapse = ";"
)
# Single OSRM /table call for all 101 x 101 pairs
osrm_url <- paste0(
"http://router.project-osrm.org/table/v1/driving/",
coords_str,
"?annotations=distance"
)
message(" Sending OSRM table request...")
## 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))
}
## OSRM unavailable — using straight-line distances as fallback.
# Replace NA / Inf with large penalty
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
# 6. MST-DFS tour
message("Solving TSP via MST-DFS...")
## 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)
## Warning: The `neimode` argument of `dfs()` is deprecated as of igraph 1.3.0.
## ℹ Please use the `mode` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
tour <- c(as.integer(dfs_res$order), 1L)
# 7. 2-opt improvement
message("Applying 2-opt improvement...")
## 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)
# 8. Stop numbering
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("\nTour length: %.1f km\n", tour_length_km))
##
## Tour length: 1978.1 km
# 9. Route geometry via OSRM /route
message("Fetching road geometry for each tour leg from OSRM...")
## 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))
}
## Leg 20 / 101 done
## Leg 40 / 101 done
## Leg 60 / 101 done
## Leg 80 / 101 done
## Leg 100 / 101 done
route_sfc <- st_sfc(leg_geoms, crs = 4326) |>
st_transform(utm32)
route_sf <- st_sf(geometry = route_sfc)
message(" Route geometry complete.")
## Route geometry complete.
# 10. Download road layer for background display
message("Downloading road background layer from Geofabrik...")
## 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
})
## Road background unavailable: The download operation was aborted. If this was not intentional, you may want to increase the timeout for internet operations to a value >= 300 by using options(timeout = ...) before re-running this function. We also suggest you to remove the partially downloaded file by running the following code (possibly in a new R session): file.remove("C:/Users/User/AppData/Local/Temp/Rtmpeo6E5S/geofabrik_denmark-latest.osm.pbf")
# 11. Plotting data frames
towns_df <- all_nodes[!all_nodes$depot,]
depot_df <- all_nodes[ all_nodes$depot,]
# 12. Build map
message("Building map...")
## 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)
)
## Warning in element_text(size = 11, face = "bold", colour = "black", hjust = 0.5, : `...` must be empty.
## ✖ Problematic argument:
## • decoration = "underline"
# 13. Show in RStudio Plots panel
print(p)
# 14. Save A4 Landscape PDF
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()
## png
## 2
cat(sprintf("\nPDF saved to:\n %s\n", pdf_file))
##
## PDF saved to:
## C:/Users/User/Desktop/CLASS NOTES/4.2 Notes/GIS Applications/Assignments/Assignment 5/denmark_tsp_route.pdf
# 15. Tour order to console
cat("\n--- Tour order ---\n")
##
## --- Tour order ---
for (i in seq_along(unique_tour)) {
idx <- unique_tour[i]
cat(sprintf(" Stop %3d : %s\n", i - 1L,
ifelse(all_nodes$depot[idx], "Copenhagen (Depot)",
all_nodes$name[idx])))
}
## Stop 0 : Copenhagen (Depot)
## Stop 1 : Town_086
## Stop 2 : Town_032
## Stop 3 : Town_006
## Stop 4 : Town_071
## Stop 5 : Town_001
## Stop 6 : Town_030
## Stop 7 : Town_011
## Stop 8 : Town_087
## Stop 9 : Town_038
## Stop 10 : Town_096
## Stop 11 : Town_066
## Stop 12 : Town_023
## Stop 13 : Town_075
## Stop 14 : Town_094
## Stop 15 : Town_031
## Stop 16 : Town_095
## Stop 17 : Town_012
## Stop 18 : Town_058
## Stop 19 : Town_039
## Stop 20 : Town_062
## Stop 21 : Town_097
## Stop 22 : Town_033
## Stop 23 : Town_015
## Stop 24 : Town_056
## Stop 25 : Town_048
## Stop 26 : Town_057
## Stop 27 : Town_089
## Stop 28 : Town_084
## Stop 29 : Town_019
## Stop 30 : Town_003
## Stop 31 : Town_065
## Stop 32 : Town_014
## Stop 33 : Town_054
## Stop 34 : Town_002
## Stop 35 : Town_085
## Stop 36 : Town_041
## Stop 37 : Town_100
## Stop 38 : Town_064
## Stop 39 : Town_008
## Stop 40 : Town_025
## Stop 41 : Town_021
## Stop 42 : Town_093
## Stop 43 : Town_074
## Stop 44 : Town_020
## Stop 45 : Town_051
## Stop 46 : Town_042
## Stop 47 : Town_079
## Stop 48 : Town_035
## Stop 49 : Town_080
## Stop 50 : Town_068
## Stop 51 : Town_007
## Stop 52 : Town_092
## Stop 53 : Town_037
## Stop 54 : Town_026
## Stop 55 : Town_076
## Stop 56 : Town_098
## Stop 57 : Town_069
## Stop 58 : Town_059
## Stop 59 : Town_088
## Stop 60 : Town_083
## Stop 61 : Town_099
## Stop 62 : Town_013
## Stop 63 : Town_091
## Stop 64 : Town_044
## Stop 65 : Town_050
## Stop 66 : Town_052
## Stop 67 : Town_004
## Stop 68 : Town_063
## Stop 69 : Town_040
## Stop 70 : Town_045
## Stop 71 : Town_018
## Stop 72 : Town_061
## Stop 73 : Town_073
## Stop 74 : Town_072
## Stop 75 : Town_081
## Stop 76 : Town_078
## Stop 77 : Town_036
## Stop 78 : Town_047
## Stop 79 : Town_090
## Stop 80 : Town_049
## Stop 81 : Town_027
## Stop 82 : Town_070
## Stop 83 : Town_043
## Stop 84 : Town_082
## Stop 85 : Town_028
## Stop 86 : Town_022
## Stop 87 : Town_005
## Stop 88 : Town_016
## Stop 89 : Town_053
## Stop 90 : Town_024
## Stop 91 : Town_029
## Stop 92 : Town_009
## Stop 93 : Town_055
## Stop 94 : Town_060
## Stop 95 : Town_010
## Stop 96 : Town_046
## Stop 97 : Town_067
## Stop 98 : Town_077
## Stop 99 : Town_017
## Stop 100 : Town_034
cat(" --> Return to Copenhagen\n")
## --> Return to Copenhagen