First we define the API key, install and load the required packages.
#Install the required packages
install.packages("rmapzen")
install.packages("leaflet")
install.packages("rgdal")
install.packages("knitr")
library(rmapzen)
library(leaflet)
library(rgdal)
library(knitr)
# To get an API key you need to visit - mapzen.com
Sys.setenv(MAPZEN_KEY = "mapzen-xxxxxx")
Next we read in a list of retail centers; which have been given a travel (minutes) threshold relative to their ranked size. Each center has an ID, name and spatial reference.
# Read the list of retail data
retail_cent <- read.csv("Comaprison_Destinations.csv")
#Plot the table
options(knitr.table.format = "html")
kable(head(retail_cent))
| cluster_id | name | rank1 | threshold | lon | lat |
|---|---|---|---|---|---|
| RC0627 | Abbotsinch Retail Park Paisley | 2 | 10 | -4.418374 | 55.86056 |
| TC0437 | Aberdeen | 5 | 30 | -2.097426 | 57.14700 |
| TC0761 | Abergavenny | 3 | 15 | -3.018699 | 51.82279 |
| TC1213 | Aberystwyth | 3 | 15 | -4.083511 | 52.41496 |
| TC0925 | Accrington | 3 | 15 | -2.364082 | 53.75337 |
| RC0226 | Aintree Shopping Park (Racecourse Retail Park) Liverpool | 3 | 15 | -2.949661 | 53.48024 |
Now we know that we know the number of retail centers, we can estimate the total cost of using Mapzen service to calculate the isochrones; which is currently set at $1.50 per 1,000 requests; under the assumption that your 100 free credits have been used up. Bargain!
# Price in $
(1.50 /1000) * nrow(retail_cent)
## [1] 1.4715
We can run a single example as follows.
#Choose a random retail centre
x <- sample(1:nrow(retail_cent),1)
# Create a lat lon pair
loc_tmp <- mz_location(retail_cent[x,"lat"],retail_cent[x,"lon"])
# Get the number of minutes
dist_tmp <- retail_cent[x,"threshold"]
# Calculate the isochrone using the Mapzen API
isos_tmp <- mz_isochrone(
loc_tmp,
costing_model = mz_costing$auto(),
contours = mz_contours(dist_tmp),
polygons = TRUE
)
And then create a map to view the isochrone.
leaflet(as_sp(isos_tmp)) %>%
addProviderTiles("CartoDB.DarkMatter") %>%
addPolygons(color = ~color, weight = 1)
Next we can loop through each of the retail centers in turn, and send the coordinate pairs to the API; pulling down each isochrone and saving this as a Shapefile.
for (x in 1:nrow(retail_cent)){
# Create a lat lon pair
loc_tmp <- mz_location(retail_cent[x,"lat"],retail_cent[x,"lon"])
# Get the number of minutes
dist_tmp <- retail_cent[x,"threshold"]
isos_tmp <- mz_isochrone(
loc_tmp,
costing_model = mz_costing$auto(),
contours = mz_contours(dist_tmp),
polygons = TRUE
)
#Convert to an SP object
isos_tmp <- as_sp(isos_tmp)
#Add ID Column; and remove all others
isos_tmp@data$cluster_id <- retail_cent[x,"cluster_id"]
isos_tmp@data <- data.frame(cluster_id=isos_tmp@data[,"cluster_id"])
# Write out the shapefile
writeOGR(isos_tmp, layer = retail_cent[x,"cluster_id"], 'Shapefiles', driver="ESRI Shapefile")
rm(list = c('loc_tmp','dist_tmp','isos_tmp'))
Sys.sleep(0.5)
}