Module 16 Spatial data and maps
Spatial data, also known as geospatial data, is a term used to describe any data related to or containing information about a specific location on a surface (often a map). We will consider distance matrix calculations for finding the shortest distance/travel time between a set of origins and destinations, geocoding which is the process of converting an address to geographic coordinates (latitude, longitude) and reverse geocoding which is the opposite process of converting a location as described by geographic coordinates (latitude, longitude) to a human-readable address or place name. Finally, we consider how to display spatial data on a map.
16.1 Learning outcomes
By the end of this module, you are expected to be able to:
- Calculate euclidean, manhattan, etc. distances.
- Calculate a distance matrix, i.e. shortest paths between places using the Google and Bing API.
- Geocode an address using the Google and Bing API.
- Add markers and routes (lines) on a map.
16.2 Services for obtatining spatial data
Often you need to connect to a service using an API for obtaining spatial data. The most common is Google and Bing (Microsoft) and to use the services you need an API key. Another service is also Here.
If you use Google Maps you can obtain an API key here. Modest to light use is free; however, you need a valid credit card. Note you must enable the APIs you intend to use. Google in fact has several services for geo-related solutions. For example, the Maps Static API provides map images, while the Geocoding API provides geocoding and reverse geocoding services. You need to enable the APIs before you use them. You will only need to do that once, and then they will be ready for you to use. Enabling the APIs just means clicking a few radio buttons on the Google Maps Platform web interface.
We will be using the ggmap package for Google services. You can add the API key using:
Note the key is stored in the environment object GGMAP_GOOGLE_API_KEY. Moreover, it is saved in the file .Renviron so that it is automatically reloaded when you restart R.
If you use Bing Maps you can obtain an API key here. No credit card is needed. You can add the API key using:
usethis::edit_r_environ() # opens the .Renviron file
Sys.setenv(BING_MAPS_API_KEY=[your key]) # so you don't have to restart RAdd the line BING_MAPS_API_KEY=[your key] and save the file.
Note your API keys is private and unique to you, so be careful not to share it online!
16.3 Calculating distances
If you need to calculate euclidean, manhattan, etc. distances, you can use the dist R function:
coord <- matrix(c(0,0, 0,1, 1,0, 1,1), ncol = 2, byrow = TRUE)
coord
#> [,1] [,2]
#> [1,] 0 0
#> [2,] 0 1
#> [3,] 1 0
#> [4,] 1 1
dist(coord)
#> 1 2 3
#> 2 1.00
#> 3 1.00 1.41
#> 4 1.41 1.00 1.00
as.matrix(dist(coord))
#> 1 2 3 4
#> 1 0.00 1.00 1.00 1.41
#> 2 1.00 0.00 1.41 1.00
#> 3 1.00 1.41 0.00 1.00
#> 4 1.41 1.00 1.00 0.00However, euclidean distances are often a poor approximation of shortest path lengths.
16.3.1 Using Google Maps
Remember to have set your API key and activate the Distance Matrix API service on the Google Cloud Platform. We consider the following places:
library(ggmap)
#> ℹ Google's Terms of Service: <https://mapsplatform.google.com>
#> Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
#> OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
#> ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
library(tidyverse)
dat <- tibble::tribble(
~Id, ~Shop, ~Address,
1L, "Bilka Esbjerg", "Stormgade 157, 6715 Esbjerg, Denmark",
2L, "Bilka Herning", "Golfvej 5, 7400 Herning, Denmark",
3L, "Bilka Hillerød", "Slotsarkaderne 26, 3400 Hillerød, Denmark",
4L, "Bilka Hjørring", "A.F Heidemannsvej 20, 9800 Hjørring, Denmark",
5L, "Bilka Holstebro", "Nyholmvej 20, 7500 Holstebro, Denmark",
)
dat
#> # A tibble: 5 × 3
#> Id Shop Address
#> <int> <chr> <chr>
#> 1 1 Bilka Esbjerg Stormgade 157, 6715 Esbjerg, Denmark
#> 2 2 Bilka Herning Golfvej 5, 7400 Herning, Denmark
#> 3 3 Bilka Hillerød Slotsarkaderne 26, 3400 Hillerød, Denmark
#> 4 4 Bilka Hjørring A.F Heidemannsvej 20, 9800 Hjørring, Denmark
#> 5 5 Bilka Holstebro Nyholmvej 20, 7500 Holstebro, DenmarkTo calculate a distance matrix we use the mapdist function:
mapdist(from = "Stormgade 157, 6715 Esbjerg, Denmark",
to = "Golfvej 5, 7400 Herning, Denmark")
#> ℹ <https://maps.googleapis.com/maps/api/distancematrix/json?origins=Stormgade+157,+6715+Esbjerg,+Denmark&destinations=Golfvej+5,+7400+Herning,+Denmark&key=xxx-RC0wa23c&mode=driving>
#> # A tibble: 1 × 9
#> from to m km miles seconds minutes hours mode
#> <chr> <chr> <int> <dbl> <dbl> <int> <dbl> <dbl> <chr>
#> 1 Stormgade 157, 6715 Esbjerg, Denmark Golfvej 5, 740… 87793 87.8 54.6 4368 72.8 1.21 driv…Note Google returns results for the fastest path between the two points. Let us try to define a function which calculate all the distances:
#' Calculate the distance matrix in long format.
#'
#' @param address A vector of addresses.
#' @param mode Driving, bicycling, walking, or transit.
#' @param symmetric Use symmetric distances (half the number of queries).
#' @return A data frame with the results
#' @note The API returns results for the fastest route.
goo_calc_distances <- function(address, mode = "driving", symmetric = TRUE) {
datDist <- expand_grid(id_from = 1:length(address), id_to = 1:length(address))
datDist <- datDist |> filter(id_from != id_to)
if (symmetric) datDist <- datDist |> filter(id_from < id_to)
datDist <- datDist |> mutate(from = address[id_from], to = address[id_to])
res <- mapdist(from = datDist |> pull(from), to = datDist |> pull(to), mode = mode)
return(left_join(datDist, res, by = c("from", "to")))
}
datDistGoo <- goo_calc_distances(dat$Address)
#> ℹ <https://maps.googleapis.com/maps/api/distancematrix/json?origins=A.F+Heidemannsvej+20,+9800+Hj%C3%B8rring,+Denmark&destinations=Nyholmvej+20,+7500+Holstebro,+Denmark&key=xxx-RC0wa23c&mode=driving>
#> ℹ <https://maps.googleapis.com/maps/api/distancematrix/json?origins=Golfvej+5,+7400+Herning,+Denmark&destinations=Slotsarkaderne+26,+3400+Hiller%C3%B8d,+Denmark%7CA.F+Heidemannsvej+20,+9800+Hj%C3%B8rring,+Denmark%7CNyholmvej+20,+7500+Holstebro,+Denmark&key=xxx-RC0wa23c&mode=driving>
#> ℹ <https://maps.googleapis.com/maps/api/distancematrix/json?origins=Slotsarkaderne+26,+3400+Hiller%C3%B8d,+Denmark&destinations=A.F+Heidemannsvej+20,+9800+Hj%C3%B8rring,+Denmark%7CNyholmvej+20,+7500+Holstebro,+Denmark&key=xxx-RC0wa23c&mode=driving>
#> ℹ <https://maps.googleapis.com/maps/api/distancematrix/json?origins=Stormgade+157,+6715+Esbjerg,+Denmark&destinations=Golfvej+5,+7400+Herning,+Denmark%7CSlotsarkaderne+26,+3400+Hiller%C3%B8d,+Denmark%7CA.F+Heidemannsvej+20,+9800+Hj%C3%B8rring,+Denmark%7CNyholmvej+20,+7500+Holstebro,+Denmark&key=xxx-RC0wa23c&mode=driving>
datDistGoo
#> # A tibble: 10 × 11
#> id_from id_to from to m km miles seconds minutes hours mode
#> <int> <int> <chr> <chr> <int> <dbl> <dbl> <int> <dbl> <dbl> <chr>
#> 1 1 2 Stormgade 157, 6715 Esbjerg, … Golf… 87793 87.8 54.6 4368 72.8 1.21 driv…
#> 2 1 3 Stormgade 157, 6715 Esbjerg, … Slot… 322365 322. 200. 11830 197. 3.29 driv…
#> 3 1 4 Stormgade 157, 6715 Esbjerg, … A.F … 322094 322. 200. 11856 198. 3.29 driv…
#> 4 1 5 Stormgade 157, 6715 Esbjerg, … Nyho… 103780 104. 64.5 5343 89.0 1.48 driv…
#> 5 2 3 Golfvej 5, 7400 Herning, Denm… Slot… 333816 334. 207. 12177 203. 3.38 driv…
#> 6 2 4 Golfvej 5, 7400 Herning, Denm… A.F … 175960 176. 109. 7934 132. 2.20 driv…
#> 7 2 5 Golfvej 5, 7400 Herning, Denm… Nyho… 43980 44.0 27.3 1798 30.0 0.499 driv…
#> 8 3 4 Slotsarkaderne 26, 3400 Hille… A.F … 377557 378. 235. 17494 292. 4.86 driv…
#> 9 3 5 Slotsarkaderne 26, 3400 Hille… Nyho… 373218 373. 232. 13802 230. 3.83 driv…
#> 10 4 5 A.F Heidemannsvej 20, 9800 Hj… Nyho… 181192 181. 113. 7996 133. 2.22 driv…Note that only the calculated distances are returned. If you want to have the whole distance matrix we define function:
#' Convert the data frame returned from calling a `___calc_distances` function to a distance matrix.
#'
#' @param dat The data frame returned from calling `___calc_distances`.
#' @param value_col The column containing the distances.
#' @return The distance matrix
as_dist_matrix <- function(dat, value_col) {
lgt <- max(dat$id_from, dat$id_to)
distanceMat<-matrix(NA, nrow=lgt, ncol = lgt)
diag(distanceMat) <- 0
map(1:nrow(dat), function(r) {
distanceMat[dat$id_from[r], dat$id_to[r]] <<- dat[[value_col]][r]
})
idx <- which(is.na(distanceMat), arr.ind = TRUE)
map(1:nrow(idx), function(r) {
distanceMat[idx[r, "row"], idx[r, "col"]] <<- distanceMat[idx[r, "col"], idx[r, "row"]]
})
return(distanceMat)
}
as_dist_matrix(datDistGoo, "km") # distances in km
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.0 87.8 322 322 104
#> [2,] 87.8 0.0 334 176 44
#> [3,] 322.4 333.8 0 378 373
#> [4,] 322.1 176.0 378 0 181
#> [5,] 103.8 44.0 373 181 0
as_dist_matrix(datDistGoo, "seconds") # travel time in seconds
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0 4368 11830 11856 5343
#> [2,] 4368 0 12177 7934 1798
#> [3,] 11830 12177 0 17494 13802
#> [4,] 11856 7934 17494 0 7996
#> [5,] 5343 1798 13802 7996 016.3.2 Using Bing Maps
Remember to have set your API key. We first define a bing_mapdist function:
library(jsonlite)
#' Distance between two locations.
#'
#' @param from From address.
#' @param to To address.
#' @param mode Walking, driving or transit.
#' @param optimize Optimize either `distance` or `time`.
#' @return
bing_mapdist <- function(from, to, mode = "driving", optimize = "time") {
if (is.data.frame(from)) {
stopifnot(all(c("from", "to") %in% names(from)))
from_to_df <- from |> select("from", "to") |> as_tibble()
}
else {
from_to_df <- tibble(from = from, to = to)
}
dat <- map_dfr(1:nrow(from_to_df), function(r) {
url <-
str_c("http://dev.virtualearth.net/REST/V1/Routes/", mode, "?wp.0=", from_to_df$from[r], "&wp.1=", from_to_df$to[r],
"&optimize=", optimize, "&ra=excludeItinerary&rpo=None&ig=false&du=km",
"&avoid=minimizeTolls&key=", Sys.getenv("BING_MAPS_API_KEY"))
url <- URLencode(url)
lst <- fromJSON(url)
if (lst$statusCode != 200) return(tibble(km = NA, seconds = NA, mode = NA))
df <- lst$resourceSets$resources[[1]]
return(tibble(km = df$travelDistance, seconds = df$travelDuration, mode = df$travelMode))
})
return(bind_cols(from_to_df, dat))
}
bing_mapdist(from = "Stormgade 157, 6715 Esbjerg, Denmark",
to = "A.F Heidemannsvej 20, 9800 Hjørring, Denmark")
#> # A tibble: 1 × 5
#> from to km seconds mode
#> <chr> <chr> <dbl> <int> <chr>
#> 1 Stormgade 157, 6715 Esbjerg, Denmark A.F Heidemannsvej 20, 9800 Hjørring, Den… 322. 10332 Driv…
bing_mapdist(from = "Stormgade 157, 6715 Esbjerg, Denmark",
to = "A.F Heidemannsvej 20, 9800 Hjørring, Denmark",
optimize = "distance")
#> # A tibble: 1 × 5
#> from to km seconds mode
#> <chr> <chr> <dbl> <int> <chr>
#> 1 Stormgade 157, 6715 Esbjerg, Denmark A.F Heidemannsvej 20, 9800 Hjørring, Den… 259. 13398 Driv…Note we here can get both the shortest and fastest path between two points. Next, we use the bing_mapdist function for calculating distances:
#' Calculate the distance matrix in long format.
#'
#' @param address A vector of addresses.
#' @param symmetric Use symmetric distances (only half the number of queries).
#' @param ... Further parameters passed to `bing_mapdist`.
#' @return A data frame with the results
#' @note The API returns results for the fastest route.
bing_calc_distances <- function(address, symmetric = TRUE, ...) {
datDist <- expand_grid(id_from = 1:length(address), id_to = 1:length(address))
datDist <- datDist |> filter(id_from != id_to)
if (symmetric) datDist <- datDist |> filter(id_from < id_to)
datDist <- datDist |> mutate(from = address[id_from], to = address[id_to])
res <- bing_mapdist(datDist |> select(from, to), ...)
return(bind_cols(datDist |> select(id_from, id_to), res))
}
datDistBing <- bing_calc_distances(dat$Address)
datDistBing
#> # A tibble: 10 × 7
#> id_from id_to from to km seconds mode
#> <int> <int> <chr> <chr> <dbl> <int> <chr>
#> 1 1 2 Stormgade 157, 6715 Esbjerg, Denmark Golfvej 5, 7400 H… 89.9 4266 Driv…
#> 2 1 3 Stormgade 157, 6715 Esbjerg, Denmark Slotsarkaderne 26… 322. 10562 Driv…
#> 3 1 4 Stormgade 157, 6715 Esbjerg, Denmark A.F Heidemannsvej… 322. 10332 Driv…
#> 4 1 5 Stormgade 157, 6715 Esbjerg, Denmark Nyholmvej 20, 750… 104. 5251 Driv…
#> 5 2 3 Golfvej 5, 7400 Herning, Denmark Slotsarkaderne 26… 334. 10822 Driv…
#> 6 2 4 Golfvej 5, 7400 Herning, Denmark A.F Heidemannsvej… 178. 7641 Driv…
#> 7 2 5 Golfvej 5, 7400 Herning, Denmark Nyholmvej 20, 750… 44.7 1856 Driv…
#> 8 3 4 Slotsarkaderne 26, 3400 Hillerød, Denmark A.F Heidemannsvej… 378. 16834 Driv…
#> 9 3 5 Slotsarkaderne 26, 3400 Hillerød, Denmark Nyholmvej 20, 750… 376. 12251 Driv…
#> 10 4 5 A.F Heidemannsvej 20, 9800 Hjørring, Denmark Nyholmvej 20, 750… 186. 7845 Driv…Note in general results from the API services are based on a set of assumptions and hence the results from different services may not be the same. You may use the different services to check if the distances are correct, by checking the results from different services:
# Difference in minutes
abs(as_dist_matrix(datDistGoo, "seconds") - as_dist_matrix(datDistBing, "seconds"))/60
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.00 1.700 21.1 25.40 1.533
#> [2,] 1.70 0.000 22.6 4.88 0.967
#> [3,] 21.13 22.583 0.0 11.00 25.850
#> [4,] 25.40 4.883 11.0 0.00 2.517
#> [5,] 1.53 0.967 25.9 2.52 0.000
# Km from 1 to 3
mapdist(from = "Stormgade 157, 6715 Esbjerg, Denmark",
to = "Slotsarkaderne 26, 3400 Hillerød, Denmark")
#> ℹ <https://maps.googleapis.com/maps/api/distancematrix/json?origins=Stormgade+157,+6715+Esbjerg,+Denmark&destinations=Slotsarkaderne+26,+3400+Hiller%C3%B8d,+Denmark&key=xxx-RC0wa23c&mode=driving>
#> # A tibble: 1 × 9
#> from to m km miles seconds minutes hours mode
#> <chr> <chr> <int> <dbl> <dbl> <int> <dbl> <dbl> <chr>
#> 1 Stormgade 157, 6715 Esbjerg, Denmark Slotsarkadern… 322365 322. 200. 11830 197. 3.29 driv…
bing_mapdist(from = "Stormgade 157, 6715 Esbjerg, Denmark",
to = "Slotsarkaderne 26, 3400 Hillerød, Denmark")
#> # A tibble: 1 × 5
#> from to km seconds mode
#> <chr> <chr> <dbl> <int> <chr>
#> 1 Stormgade 157, 6715 Esbjerg, Denmark Slotsarkaderne 26, 3400 Hillerød, Denmark 322. 10562 Driv…16.4 Geocoding and reverse geocoding
16.4.1 Using Google maps
Remember to enable the Geocoding API. To get the coordinates of an address we use the geocode function:
geocode("Stormgade 157, 6715 Esbjerg, Denmark")
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Stormgade+157,+6715+Esbjerg,+Denmark&key=xxx-RC0wa23c>
#> # A tibble: 1 × 2
#> lon lat
#> <dbl> <dbl>
#> 1 8.46 55.5We may use mutate_geocode to add the coordinates to a dataset:
dat <- dat |>
mutate_geocode(Address)
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Stormgade+157,+6715+Esbjerg,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Golfvej+5,+7400+Herning,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Slotsarkaderne+26,+3400+Hiller%C3%B8d,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=A.F+Heidemannsvej+20,+9800+Hj%C3%B8rring,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Nyholmvej+20,+7500+Holstebro,+Denmark&key=xxx-RC0wa23c>
dat
#> # A tibble: 5 × 5
#> Id Shop Address lon lat
#> <int> <chr> <chr> <dbl> <dbl>
#> 1 1 Bilka Esbjerg Stormgade 157, 6715 Esbjerg, Denmark 8.46 55.5
#> 2 2 Bilka Herning Golfvej 5, 7400 Herning, Denmark 9.01 56.1
#> 3 3 Bilka Hillerød Slotsarkaderne 26, 3400 Hillerød, Denmark 12.3 55.9
#> 4 4 Bilka Hjørring A.F Heidemannsvej 20, 9800 Hjørring, Denmark 10.0 57.4
#> 5 5 Bilka Holstebro Nyholmvej 20, 7500 Holstebro, Denmark 8.62 56.4To reverse geocode use the revgeocode function:
revgeocode(c(dat$lon[1], dat$lat[1]))
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?latlng=55.5107577,8.4559613&key=xxx-RC0wa23c>
#> Warning: Multiple addresses found, the first will be returned:
#> ! Stormgade 157, 6715 Esbjerg Kommune, Denmark
#> ! Stormgade 157, 6715 Esbjerg N, Denmark
#> ! GF64+89 Esbjerg, Denmark
#> ! Stormgade 156, 6715 Esbjerg Kommune, Denmark
#> ! Stormgade 158-155, 6715 Esbjerg Kommune, Denmark
#> ! 6715 Esbjerg N, Denmark
#> ! Esbjerg N, 6715 Esbjerg N, Denmark
#> ! Esbjerg N, 6715 Esbjerg, Denmark
#> ! Esbjerg, Denmark
#> ! Esbjerg Municipality, Denmark
#> ! Region of Southern Denmark, Denmark
#> ! Denmark
#> [1] "Stormgade 157, 6715 Esbjerg Kommune, Denmark"To apply it to the dataset use:
dat <- dat |>
mutate(AddressGeoGoo =
map_chr(1:n(), function(i) {
revgeocode(c(lon = lon[i], lat = lat[i]))
}))
#> ! Stormgade 157, 6715 Esbjerg Kommune, Denmark
#> ! Stormgade 157, 6715 Esbjerg N, Denmark
#> ! GF64+89 Esbjerg, Denmark
#> ! Stormgade 156, 6715 Esbjerg Kommune, Denmark
#> ! Stormgade 158-155, 6715 Esbjerg Kommune, Denmark
#> ! 6715 Esbjerg N, Denmark
#> ! Esbjerg N, 6715 Esbjerg N, Denmark
#> ! Esbjerg N, 6715 Esbjerg, Denmark
#> ! Esbjerg, Denmark
#> ! Esbjerg Municipality, Denmark
#> ! Region of Southern Denmark, Denmark
#> ! Denmark
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?latlng=56.1357174,9.0053714&key=xxx-RC0wa23c>
#> ! Golfvej 5, 7400 Herning, Denmark
#> ! 42P4+74 Herning, Denmark
#> ! Engdahlsvej, 7400 Herning, Denmark
#> ! 7400 Herning, Denmark
#> ! Herning, 7400 Herning, Denmark
#> ! Herning Municipality, Denmark
#> ! Central Denmark Region, Denmark
#> ! Denmark
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?latlng=55.9288332,12.2978145&key=xxx-RC0wa23c>
#> ! Slotsarkaderne 30, 3400 Hillerød, Denmark
#> ! Slotsarkaderne 26, 3400 Hillerød, Denmark
#> ! W7HX+G4 Hillerød, Denmark
#> ! Hillerød, 3400 Hillerød, Denmark
#> ! 3400 Hillerød, Denmark
#> ! Hillerød Municipality, Denmark
#> ! Capital Region of Denmark, Denmark
#> ! Denmark
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?latlng=57.4441555,10.0034928&key=xxx-RC0wa23c>
#> ! A F Heidemanns Vej 20, 9800 Hjørring, Denmark
#> ! A F Heidemanns Vej 21, 9800 Hjørring, Denmark
#> ! C2V3+M9 Hjørring, Denmark
#> ! 9800 Hjørring, Denmark
#> ! Hjørring, 9800 Hjørring, Denmark
#> ! Hjørring Municipality, Denmark
#> ! North Denmark Region, Denmark
#> ! Denmark
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?latlng=56.3779116,8.6222409&key=xxx-RC0wa23c>
#> ! Nyholmvej 20, 7500 Holstebro, Denmark
#> ! Nyholmvej 20c, 7500 Holstebro, Denmark
#> ! 9JHC+5V Holstebro, Denmark
#> ! Nyholmvej 15, 7500 Holstebro, Denmark
#> ! Nyholmvej 19-17, 7500 Holstebro, Denmark
#> ! 7500 Holstebro, Denmark
#> ! Holstebro, 7500 Holstebro, Denmark
#> ! Holstebro Municipality, Denmark
#> ! Central Denmark Region, Denmark
#> ! Denmark
#> Warning: There were 5 warnings in `mutate()`.
#> The first warning was:
#> ℹ In argument: `AddressGeoGoo = map_chr(...)`.
#> Caused by warning:
#> ! Multiple addresses found, the first will be returned:
#> ℹ Run `dplyr::last_dplyr_warnings()` to see the 4 remaining warnings.
dat
#> # A tibble: 5 × 6
#> Id Shop Address lon lat AddressGeoGoo
#> <int> <chr> <chr> <dbl> <dbl> <chr>
#> 1 1 Bilka Esbjerg Stormgade 157, 6715 Esbjerg, Denmark 8.46 55.5 Stormgade 157, 671…
#> 2 2 Bilka Herning Golfvej 5, 7400 Herning, Denmark 9.01 56.1 Golfvej 5, 7400 He…
#> 3 3 Bilka Hillerød Slotsarkaderne 26, 3400 Hillerød, Denmark 12.3 55.9 Slotsarkaderne 30,…
#> 4 4 Bilka Hjørring A.F Heidemannsvej 20, 9800 Hjørring, Denmark 10.0 57.4 A F Heidemanns Vej…
#> 5 5 Bilka Holstebro Nyholmvej 20, 7500 Holstebro, Denmark 8.62 56.4 Nyholmvej 20, 7500…16.4.2 Use Bing Maps
To get the coordinates of an address we define the bing_geocode function:
#' Geocode addresses
#'
#' @param address Address(es) to geocode.
#' @return The coordinates as a data frame.
bing_geocode <- function(address) {
dat <- map_dfr(address, function(s) {
url <- str_c("http://dev.virtualearth.net/REST/v1/Locations?q=",
s, "&key=", Sys.getenv("BING_MAPS_API_KEY"))
url <- URLencode(url)
lst <- fromJSON(url)
if (lst$statusCode != 200) return(tibble(lon = NA, lat = NA))
v <- lst$resourceSets$resources[[1]]$point$coordinates[[1]]
if (is.null(v)) return(tibble(lon = NA, lat = NA))
return(tibble(lon = v[2], lat = v[1]))
})
return(dat)
}
bing_geocode("A.F. Heidemannsvej 20, 9800 Hjørring, Denmark")
#> # A tibble: 1 × 2
#> lon lat
#> <dbl> <dbl>
#> 1 10.0 57.4
bing_geocode(dat$Address)
#> # A tibble: 5 × 2
#> lon lat
#> <dbl> <dbl>
#> 1 8.46 55.5
#> 2 9.00 56.1
#> 3 12.3 55.9
#> 4 10.0 57.4
#> 5 8.62 56.4We may use bing_mutate_geocode to add the coordinates to a dataset:
dat <- dat |> select(-lat, -lon)
mutate_bing_geocode <- function (data, address, ...) {
adr <- data[[deparse(substitute(Address))]]
gcdf <- bing_geocode(adr)
bind_cols(data, gcdf)
}
dat <- dat |>
mutate_bing_geocode(Address)
dat
#> # A tibble: 5 × 6
#> Id Shop Address AddressGeoGoo lon lat
#> <int> <chr> <chr> <chr> <dbl> <dbl>
#> 1 1 Bilka Esbjerg Stormgade 157, 6715 Esbjerg, Denmark Stormgade 157, 671… 8.46 55.5
#> 2 2 Bilka Herning Golfvej 5, 7400 Herning, Denmark Golfvej 5, 7400 He… 9.00 56.1
#> 3 3 Bilka Hillerød Slotsarkaderne 26, 3400 Hillerød, Denmark Slotsarkaderne 30,… 12.3 55.9
#> 4 4 Bilka Hjørring A.F Heidemannsvej 20, 9800 Hjørring, Denmark A F Heidemanns Vej… 10.0 57.4
#> 5 5 Bilka Holstebro Nyholmvej 20, 7500 Holstebro, Denmark Nyholmvej 20, 7500… 8.62 56.4To reverse geocode use the bing_revgeocode function:
#' Reverse geocode coordinates
#'
#' @param coordinates A vector with two elements (lon, lat).
#' @return The address as a data frame.
bing_revgeocode <- function(coordinates) {
url <- str_c("http://dev.virtualearth.net/REST/v1/Locations/",
coordinates[2], ",", coordinates[1],
"?key=", Sys.getenv("BING_MAPS_API_KEY"))
url <- URLencode(url)
lst <- fromJSON(url)
if (lst$statusCode != 200) return(NA)
v <- lst$resourceSets$resources[[1]]$address$formattedAddress
if (is.null(v)) return(NA)
return(v)
}
bing_revgeocode(c(dat$lon[1], dat$lat[1]))
#> [1] "Stormgade 157, 6715 Esbjerg, Denmark"To apply it to the dataset use:
dat <- dat |>
mutate(AddressGeoBing =
map_chr(1:n(), function(i) {
bing_revgeocode(c(lon = lon[i], lat = lat[i]))
}))
dat
#> # A tibble: 5 × 7
#> Id Shop Address AddressGeoGoo lon lat AddressGeoBing
#> <int> <chr> <chr> <chr> <dbl> <dbl> <chr>
#> 1 1 Bilka Esbjerg Stormgade 157, 6715 Esbjerg, Denma… Stormgade 15… 8.46 55.5 Stormgade 157…
#> 2 2 Bilka Herning Golfvej 5, 7400 Herning, Denmark Golfvej 5, 7… 9.00 56.1 Golfvej 5, 74…
#> 3 3 Bilka Hillerød Slotsarkaderne 26, 3400 Hillerød, … Slotsarkader… 12.3 55.9 Slotsarkadern…
#> 4 4 Bilka Hjørring A.F Heidemannsvej 20, 9800 Hjørrin… A F Heideman… 10.0 57.4 A F Heidemann…
#> 5 5 Bilka Holstebro Nyholmvej 20, 7500 Holstebro, Denm… Nyholmvej 20… 8.62 56.4 Nyholmvej 20,…16.5 Adding markers and routes to a map
Leaflet is an open-source JavaScript library for interactive maps. It’s used by websites ranging from The New York Times and The Washington Post to GitHub and Flickr, as well as GIS specialists like OpenStreetMap, Mapbox, and CartoDB. The leaflet package makes it easy to create Leaflet maps from R.
Note Leaflet is open-source and free so you do not need an API key for making maps. If you would like to use Google maps instead then have a look at the googleway package instead.
First let us create a map with two base layers
library(leaflet)
library(htmlwidgets)
m <- leaflet() |>
# Base maps
addTiles(group = "Map") |>
addProviderTiles('Esri.WorldImagery', group = "Satelite") |>
addProviderTiles("CartoDB.PositronOnlyLabels", group = "Map") |>
# Center and zoom
setView(10.2, 56.2, zoom = 7) |>
# Layer control
addLayersControl(
baseGroups = c("Map", "Satelite"),
options = layersControlOptions(collapsed = TRUE)
)
mNext, let us add the places:
dat <- tibble::tribble(
~Id, ~Shop, ~Address,
1L, "Bilka Esbjerg", "Stormgade 157, 6715 Esbjerg, Denmark",
2L, "Bilka Herning", "Golfvej 5, 7400 Herning, Denmark",
3L, "Bilka Hillerød", "Slotsarkaderne 26, 3400 Hillerød, Denmark",
4L, "Bilka Hjørring", "A.F Heidemannsvej 20, 9800 Hjørring, Denmark",
5L, "Bilka Holstebro", "Nyholmvej 20, 7500 Holstebro, Denmark",
)
dat <- dat |> mutate_geocode(Address)
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Stormgade+157,+6715+Esbjerg,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Golfvej+5,+7400+Herning,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Slotsarkaderne+26,+3400+Hiller%C3%B8d,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=A.F+Heidemannsvej+20,+9800+Hj%C3%B8rring,+Denmark&key=xxx-RC0wa23c>
#> ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Nyholmvej+20,+7500+Holstebro,+Denmark&key=xxx-RC0wa23c>
m <- m |>
addMarkers(~lon, ~lat, popup = ~Address, label = ~str_c(Id, " - ", Shop), data = dat)
mTo add a line between a set of points use the addPolylines function:
routeIds <- c(2, 5, 1, 2)
route <- dat[routeIds,]
m |> addPolylines(lng = ~lon, lat = ~lat, data = route, weight = 2, label = "Route 1")A more advanced setup is to use the addFlows function from the leaflet.minicharts package. First let us define some lines:
datLines <- tibble(FromId = c(2, 5, 1, 2, 4, 3), ToId = c(5, 1, 2, 4, 3, 2))
datLines <- left_join(datLines, dat |> select(-Shop, -Address), by = c("FromId" = "Id")) |>
rename("FromLat" = "lat", "FromLon" = "lon") |>
left_join(dat |> select(-Shop, -Address), by = c("ToId" = "Id")) |>
rename("ToLat" = "lat", "ToLon" = "lon")
datLines
#> # A tibble: 6 × 6
#> FromId ToId FromLon FromLat ToLon ToLat
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2 5 9.01 56.1 8.62 56.4
#> 2 5 1 8.62 56.4 8.46 55.5
#> 3 1 2 8.46 55.5 9.01 56.1
#> 4 2 4 9.01 56.1 10.0 57.4
#> 5 4 3 10.0 57.4 12.3 55.9
#> 6 3 2 12.3 55.9 9.01 56.1Note we have a from/to pair in each row. We add the lines to the map:
library(leaflet.minicharts)
m |> addFlows(datLines$FromLon, datLines$FromLat, datLines$ToLon, datLines$ToLat,
flow = 1, maxFlow = 20, opacity = 0.5, popup = popupArgs(noPopup = TRUE))Let us define a function for adding routes:
add_route <- function(dat = NULL, route, solution = 1) {
route_id = if_else(is.null(dat), 1, max(dat$RouteId) + 1)
tmp <- tibble(From = route[1:(length(route)-1)], To = route[2:length(route)])
tmp <- tmp |> mutate(Sol = solution, RouteId = route_id)
dat <- bind_rows(dat, tmp)
}
datLines <- add_route(route = c(2, 5, 1, 2)) |>
add_route(route = c(2, 4, 3, 2)) |>
add_route(route = c(2, 3, 1, 2), solution = 2) |>
add_route(route = c(2, 5, 4, 2), solution = 2)
#> Warning in max(dat$RouteId): no non-missing arguments to max; returning -Inf
datLines <- left_join(datLines, dat |> select(-Shop, -Address), by = c("From" = "Id")) |>
rename("FromLat" = "lat", "FromLon" = "lon") |>
left_join(dat |> select(-Shop, -Address), by = c("To" = "Id")) |>
rename("ToLat" = "lat", "ToLon" = "lon") |>
mutate(group = str_c("Solution ", Sol))
datLines
#> # A tibble: 12 × 9
#> From To Sol RouteId FromLon FromLat ToLon ToLat group
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 2 5 1 1 9.01 56.1 8.62 56.4 Solution 1
#> 2 5 1 1 1 8.62 56.4 8.46 55.5 Solution 1
#> 3 1 2 1 1 8.46 55.5 9.01 56.1 Solution 1
#> 4 2 4 1 2 9.01 56.1 10.0 57.4 Solution 1
#> 5 4 3 1 2 10.0 57.4 12.3 55.9 Solution 1
#> 6 3 2 1 2 12.3 55.9 9.01 56.1 Solution 1
#> 7 2 3 2 3 9.01 56.1 12.3 55.9 Solution 2
#> 8 3 1 2 3 12.3 55.9 8.46 55.5 Solution 2
#> 9 1 2 2 3 8.46 55.5 9.01 56.1 Solution 2
#> 10 2 5 2 4 9.01 56.1 8.62 56.4 Solution 2
#> 11 5 4 2 4 8.62 56.4 10.0 57.4 Solution 2
#> 12 4 2 2 4 10.0 57.4 9.01 56.1 Solution 2A map with Solution 1:
col_pal <- rainbow(max(datLines$RouteId))
datP <- datLines |> filter(Sol == 1)
m |> addFlows(datP$FromLon, datP$FromLat, datP$ToLon, datP$ToLat,
flow = datP$RouteId, maxThickness = 2, color = col_pal[datP$RouteId],
opacity = 0.5, popup = popupArgs(labels = "Route"))A map with Solution 2:
col_pal <- rainbow(max(datLines$RouteId))
datP <- datLines |> filter(Sol == 2)
m |> addFlows(datP$FromLon, datP$FromLat, datP$ToLon, datP$ToLat,
flow = datP$RouteId, maxThickness = 2, color = col_pal[datP$RouteId],
opacity = 0.5, popup = popupArgs(labels = "Route"))An interactive map with both solutions:
mm <- m
res <- map(1:nrow(datLines), function(i) {
tmp <- tibble(lat = c(datLines$FromLat[i], datLines$ToLat[i]), lon = c(datLines$FromLon[i], datLines$ToLon[i]))
mm <<- mm |> addPolylines(lng = ~lon, lat = ~lat,
color = col_pal[datLines$Sol[i]], group = datLines$group[i],
data = tmp, weight = 2, label = datLines$group[i])
return(invisible(NULL))
})
mm <- mm |>
# Layer control
addLayersControl(
baseGroups = c("Map", "Satelite"),
overlayGroups = unique(datLines$group),
options = layersControlOptions(collapsed = FALSE)
) |>
hideGroup(unique(datLines$group)) |>
showGroup("Solution 1")
mmNote unfortunately arrows cannot be shown in this case.