Inspired by David Schoch’s blog post,
Traveling Beerdrinker Problem.
Check out his blog, he has some amazing posts!
Introduction
Luxembourg, as any proper European country, is full of castles. According to Wikipedia,
“By some optimistic estimates, there are as many as 130 castles in Luxembourg but more realistically
there are probably just over a hundred, although many of these could be considered large residences
or manor houses rather than castles”.
I see the editors are probably German or French, calling our castles manor houses! They only say
that because Luxembourg is small, so our castles must be small too, right?
Banter aside, with that many castles, what is the best way to visit them all? And by best way I mean
shortest way. This is a classical Travelling salesman problem. To solve this, I need the following elements:
A list of castles to visit, with their coordinates
The distances between these castles to each other
A program to solve the TSP
Let’s start by loading some packages:
library(tidyverse)
library(magrittr)
library(rvest)
library(curl)
library(brotools)
library(RJSONIO)
library(TSP)
library(ggimage)
First step; scrape the data.
Scraping the data (that’s the data science part)
Let’s start by having a list of castles. For this, I go to the French Wikipedia page of
Luxembourguish castles.
The Luxembourguish page is more exhaustive,
but the names are in Luxembourguish, and I doubt that
OpenStreetMap, which I’ll use to get the coordinates, understands Luxembourguish.
This list has around 50 castles, a reasonable amount of castles. Scraping the table is quite easy:
page %
html_table(fill = TRUE) %__%
select(Nom, Localité) %__%
mutate(query = paste0(Nom, &quot;, &quot;, Localité))
I also add a query column which concatenates the name of the castle (“Nom”) to where it is found
(“Localité”). The query should be a better choice that simply the castle name to get the coordinates.
Now, I need to add the coordinates to this data frame. For this, I use a function I found online
that gets the coordinates from OpenStreetMap:
## geocoding function using OSM Nominatim API
## details: http://wiki.openstreetmap.org/wiki/Nominatim
## made by: D.Kisler
#https://datascienceplus.com/osm-nominatim-with-r-getting-locations-geo-coordinates-by-its-address/
nominatim_osm %
mutate(geolocation = ifelse(Nom == &quot;Château de Septfontaines&quot;, &quot;5.9617443, 49.7006292&quot;, geolocation)) %__%
mutate(geolocation = ifelse(Nom == &quot;Château de Senningen&quot;, &quot;6.2342581, 49.6464632&quot;, geolocation)) %__%
mutate(geolocation = ifelse(Nom == &quot;Château de Schauwenburg&quot;, &quot;6.0478341, 49.6110245&quot;, geolocation)) %__%
mutate(geolocation = ifelse(Nom == &quot;Château de Schuttbourg&quot;, &quot;5.8980951, 49.7878706&quot;, geolocation)) %__%
mutate(geolocation = ifelse(Nom == &quot;Château de Meysembourg&quot;, &quot;6.1864882, 49.7704348&quot;, geolocation)) %__%
mutate(geolocation = ifelse(Nom == &quot;Château de Mamer&quot;, &quot;6.0232432, 49.6262397&quot;, geolocation)) %__%
mutate(geolocation = ifelse(Nom == &quot;Château de Born&quot;, &quot;6.5125214, 49.7611168&quot;, geolocation)) %__%
# Found chateau de Betzdorf in Germany, not Luxembourg:
mutate(geolocation = ifelse(Nom == &quot;Château Betzdorf&quot;, &quot;6.330278, 49.694167&quot;, geolocation)) %__%
# Found château de Clemency in France, not Luxembourg:
mutate(geolocation = ifelse(Nom == &quot;Château de Clemency&quot;, &quot;5.874167, 49.598056&quot;, geolocation)) %__%
separate(geolocation, into = c(&quot;lon&quot;, &quot;lat&quot;), sep = &quot;,&quot;) %__%
filter(!is.na(lat)) %__%
mutate(lon = map(lon, extract_numbers)) %__%
mutate(lat = map(lat, extract_numbers)) %__%
# Château de Beaufort found is in southern France, not the one in lux
# Château de Dudelange is wrong (same as Bettembourg)
# Château de Pétange is wrong (same as Differdange)
# Château d'Urspelt is wrong (same as Clervaux)
# Château d'Hesperange is wrong (same as Palais Grand-Ducal)
mutate(lon = ifelse(Nom == &quot;Château de Beaufort&quot;, &quot;6.2865176&quot;, lon),
lat = ifelse(Nom == &quot;Château de Beaufort&quot;, &quot;49.8335306&quot;, lat)) %__%
mutate(lon = ifelse(Nom == &quot;Château Dudelange&quot;, &quot;6.0578438&quot;, lon),
lat = ifelse(Nom == &quot;Château Dudelange&quot;, &quot;49.4905049&quot;, lat)) %__%
mutate(lon = ifelse(Nom == &quot;Château de Pétange&quot;, &quot;6.105703&quot;, lon),
lat = ifelse(Nom == &quot;Château de Pétange&quot;, &quot;49.7704746&quot;, lat)) %__%
mutate(lon = ifelse(Nom == &quot;Château d' Urspelt&quot;, &quot;6.043375&quot;, lon),
lat = ifelse(Nom == &quot;Château d' Urspelt&quot;, &quot;50.075342&quot;, lat)) %__%
mutate(lon = ifelse(Nom == &quot;Château d'Hesperange&quot;, &quot;6.1524302&quot;, lon),
lat = ifelse(Nom == &quot;Château d'Hesperange&quot;, &quot;49.573071&quot;, lat)) %__%
mutate(latlon = paste0(lat, &quot;,&quot;, lon)) %__%
mutate(lon = as.numeric(lon), lat = as.numeric(lat))
In the end, I have 48 castles, 2 of them were not found neither by OpenStreetMap nor
Google Maps.
Now I can get the distances matrix. For this, I opened an account at Graphhopper
and used their Matrix API. When you open
a free account, you get a standard account for free for two weeks, which was perfect for this little
exercise.
To use the Matrix API you can make a call with curl from your terminal, like this:
curl &quot;https://graphhopper.com/api/1/matrix?point=49.932707,11.588051&amp;point=50.241935,10.747375&amp;point=50.118817,11.983337&amp;type=json&amp;vehicle=car&amp;debug=true&amp;out_array=weights&amp;out_array=times&amp;out_array=distances&amp;key=[YOUR_KEY]&quot;
To use this from R, I use the {curl} package and the curl_download() function to download and
write the output to disk.
I built the url like this. First, the “points” part:
points