Geospatial Analysis Tutorial 2: Identifying locations for Opening Restaurents
2024-04-23
Introduction
In the first tutorial, we looked at how we can download the geojson file from the internet and then use it to visualize maps of various countries(we looked at Indonesia). In this blog, we will take a step further and identify potential location for opening an Indian restaurant for Singapore.
Step 0: What will be the steps to achieve this
- Get the geojson map for Singapore
- Get the metric at a planning area level or unit of analysis
- Shortlist the area of interest based on heat map
- Compare key locations in terms of constant drive time distances
Step 1: Importing the libraries
package.name<-c("tidyverse","sf","raster","tmap",
"leaflet","mapview","viridis","ggthemes",
"htmltools","osrm","mapboxapi")
for(i in package.name){
if(!require(i,character.only = T)){
install.packages(i)
}
library(i,character.only = T)
}Loading required package: tidyverse── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.0 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Loading required package: sf
Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
Loading required package: raster
Loading required package: sp
Attaching package: 'raster'
The following object is masked from 'package:dplyr':
select
Loading required package: tmap
Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
Loading required package: leaflet
Loading required package: mapview
Loading required package: viridis
Loading required package: viridisLite
Loading required package: ggthemes
Loading required package: htmltools
Loading required package: osrm
Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright
Routing: OSRM - http://project-osrm.org/
Loading required package: mapboxapi
Usage of the Mapbox APIs is governed by the Mapbox Terms of Service.
Please visit https://www.mapbox.com/legal/tos/ for more information.Step 2:Downloading the geo-json map for Singapore
I found that the map shared in a Kaggle challenge seems to fulfil all the requirements. You can download it based on the below steps
go to the link and then download the file as shown below
Step 3: Reading the geojson file into R
Now we can use the downloaded file and read it using read_sf function
file_json<-"district_and_planning_area.geojson"
sg_map<-read_sf(file_json)Plotting the Map to check its features
leaflet(sg_map)%>%
addTiles()%>%
addPolygons(fillColor = "white",
color = "#000000",
weight = 5,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=2,
fillOpacity = 3)We can note here is that the map has certain geographical features such as lat long details, some relief information,etc.We would need to have a vanilla version of the Map so that it acts as a whiteboard for our analysis.
Step 4: Cleaning the Map
Removing/Commenting addTiles converts the map into a plain format
leaflet(sg_map)%>%
# addTiles()%>%
addPolygons(fillColor = "white",
color = "#000000",
weight = 2,
fillOpacity = 2,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=2,
fillOpacity = 3)For getting the white background, we will be using the htmltools features
backg<-htmltools::tags$style(".leaflet-container{ background:white;}")
leaflet(sg_map)%>%
htmlwidgets::prependContent(backg)%>%
# addTiles()%>%
addPolygons(fillColor = "white",
color = "black",
weight = 2,
fillOpacity = 2,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=2,
fillOpacity = 3)Step 5: Adding some dummy data to the map
Lets add a column called Metric to the sg_map table. This column can be population, birth rate, total residential blocks etc.
sg_map2<-sg_map%>%
mutate(Metric=rnorm(nrow(sg_map)))
head(sg_map2)Simple feature collection with 6 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.7365 ymin: 1.254808 xmax: 103.9636 ymax: 1.416214
Geodetic CRS: WGS 84
# A tibble: 6 × 4
district planning_area geometry Metric
<chr> <chr> <MULTIPOLYGON [°]> <dbl>
1 East Bedok (((103.9321 1.305548, 103.9321 1.305… 1.12
2 Central Bukit Timah (((103.7977 1.348128, 103.7981 1.347… -1.17
3 West Bukit Batok (((103.7641 1.370011, 103.7644 1.369… -0.0489
4 Central Bukit Merah (((103.8236 1.260178, 103.8236 1.260… -0.704
5 North Central Water Catchment (((103.807 1.411259, 103.8075 1.4098… 0.681
6 Central Downtown Core (((103.8666 1.303862, 103.867 1.3033… 0.130 Step 6: Plotting the dummy data on the map
backg<-htmltools::tags$style(".leaflet-container{ background:white;}")
leaflet(sg_map)%>%
htmlwidgets::prependContent(backg)%>%
# addTiles()%>%
addPolygons(fillColor = "white",
color = "black",
weight = 2,
fillOpacity = 2,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolygons(data=sg_map2%>%
filter(Metric >= 0),
fillColor = "#CC0066",
# fill=T,color="Rose pink",
weight = 2,
fillOpacity = 2,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"),
labelOptions = labelOptions(noHide = T))%>%
addPolylines(color="black",
weight=2,
fillOpacity = 3)%>%
addPolygons(data=sg_map2%>%
filter(Metric < 0),
fillColor = "blue",
# fill=T,color="Rose pink",
weight = 2,
fillOpacity = 0.4,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"),
labelOptions = labelOptions(noHide = T))%>%
addPolylines(color="black",
weight=2,
fillOpacity = 3)Normally I have seen people who have worked in Tableau complain about lack of color options in R. A quick tip for this is to use powerpoint option and identify the hex code of the color using the below screenshot.
backg<-htmltools::tags$style(".leaflet-container{ background:white;}")
leaflet(sg_map)%>%
htmlwidgets::prependContent(backg)%>%
# addTiles()%>%
addPolygons(fillColor = "white",
color = "black",
weight = 2,
fillOpacity = 2,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolygons(data=sg_map2%>%
filter(Metric >= 0),
fillColor = "#CC0066",
# fill=T,color="Rose pink",
weight = 2,
fillOpacity = 2,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"),
labelOptions = labelOptions(noHide = T))%>%
addPolylines(color="black",
weight=2,
fillOpacity = 3)%>%
addPolygons(data=sg_map2%>%
filter(Metric < 0),
fillColor = "#00B0F0",
# fill=T,color="Rose pink",
weight = 2,
fillOpacity = 2,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"),
labelOptions = labelOptions(noHide = T))%>%
addPolylines(color="black",
weight=2,
fillOpacity = 3)You can then play around with the fillOpacity parameter to adjust the intensity of the color.
Step 7: Zooming into Serangoon and adjacent areas for illustration
Out of the all the planning area, lets say we want to focus on Serangoon.Lets plot Serangoon, Hougang, Sengkang and Ang Mo Kio portion of the Map only.This is done so that we can plot adjacent areas and show some references in the Map.
leaflet(sg_map%>%
filter(planning_area %in% c("Serangoon","Hougang",
"Sengkang","Ang Mo Kio")))%>%
addProviderTiles("CartoDB.Positron",group = "Greyscale")%>%
addPolygons(fillColor = "white",
color = "grey",
weight = 5,
fillOpacity = 0.3,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=5,
fillOpacity = 5)Now lets find some competing location in this region and get their 10 mins drive time coverage. Drive time coverage area is also know as Isochrone. These represents the areas reachable within a given time span from a point. These areas of equal travel time are called isochrones.
Step 8: Identifying coverage of few competing locations
The following three location can be used for illustration purposes
- NEX Mall Serangoon (103.8717,1.3516)
- Greenwich V (103.8695,1.3891)
- Kebun Baru Mall (103.8392,1.3696)
loc_list<-list(c(103.8717,1.3516),
c(103.8695,1.3891),
c(103.8392,1.3696))
loc<-do.call(rbind.data.frame,loc_list)
colnames(loc)<-c("long","lat")
loc long lat
1 103.8717 1.3516
2 103.8695 1.3891
3 103.8392 1.3696loc2<-loc%>%
mutate(Location = c("NEX Mall Serangoon",
"Greenwich V",
"Kebun Baru Mall"))%>%
mutate(Index=1:n())
loc2 long lat Location Index
1 103.8717 1.3516 NEX Mall Serangoon 1
2 103.8695 1.3891 Greenwich V 2
3 103.8392 1.3696 Kebun Baru Mall 3leaflet(sg_map%>%
filter(planning_area %in% c("Serangoon","Hougang",
"Sengkang","Ang Mo Kio")))%>%
addProviderTiles("CartoDB.Positron",group = "Greyscale")%>%
addPolygons(fillColor = "white",
color = "grey",
weight = 5,
fillOpacity = 0.3,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=5,
fillOpacity = 5)%>%
addMarkers(lng=loc2$long,lat=loc2$lat,
label=~lapply(paste0(loc2[["Location"]],sep="<br>"),HTML),
labelOptions = labelOptions(noHide=T))Step 9: Getting The catchment areas of these stores
Now since we have some competing locations, we need to compare them in terms of their catchment. One way to do that is to identify the area of influence around these locations.We can use isochrone- which is area of constant travel time around a point.In R, we can do it using mapboxapi package
For using mapboxapi package, you will have to create an account and enter card details. The link is https://walker-data.com/mapboxapi/
We have to create the location value which will be combination of long and lat and pass it to mb_isochrone function
isoc_item<-lapply(loc2[["Index"]], function(i){
df2<-loc2%>%
filter(Index==i)
nm<-c(as.character(df2$long),as.character(df2$lat))
# print(nm)
drive_10min <- mb_isochrone(nm,
profile = "driving-traffic",
time = 10)
return(drive_10min)
})Lets combine the info into a data frames
interim.df<-do.call(rbind.data.frame,isoc_item)%>%
mutate(Location=loc2$Location)
interim.dfSimple feature collection with 3 features and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 103.8022 ymin: 1.309494 xmax: 103.9199 ymax: 1.430516
Geodetic CRS: WGS 84
# A tibble: 3 × 3
time geometry Location
* <int> <POLYGON [°]> <chr>
1 10 ((103.8557 1.400749, 103.8555 1.397835, 103.8527 1.397717, 103… NEX Mal…
2 10 ((103.8735 1.430516, 103.8722 1.430383, 103.8668 1.424824, 103… Greenwi…
3 10 ((103.8102 1.40541, 103.8101 1.403541, 103.8109 1.403261, 103.… Kebun B…Step 10: Plotting the catchment areas
Creating the color palette for visualization of Isochrones.
viz.pal<-rev(topo.colors(3))
viz.pal[1] "#FFFF00" "#00FF4D" "#4C00FF"NEX Mall Serangoon
leaflet(sg_map%>%
filter(planning_area %in% c("Serangoon","Hougang",
"Sengkang","Ang Mo Kio")))%>%
addProviderTiles("CartoDB.Positron",group = "Greyscale")%>%
addPolygons(fillColor = "white",
color = "grey",
weight = 5,
fillOpacity = 0.3,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=5,
fillOpacity = 5)%>%
addMarkers(lng=loc2$long,lat=loc2$lat,
label=~lapply(paste0(loc2[["Location"]],sep="<br>"),HTML),
labelOptions = labelOptions(noHide=T))%>%
addPolygons(data=interim.df%>%
filter(Location=="NEX Mall Serangoon"),
fillColor = viz.pal[1],
color = "#000000",
weight = 4,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{Location}"))Greenwich V
leaflet(sg_map%>%
filter(planning_area %in% c("Serangoon","Hougang",
"Sengkang","Ang Mo Kio")))%>%
addProviderTiles("CartoDB.Positron",group = "Greyscale")%>%
addPolygons(fillColor = "white",
color = "grey",
weight = 5,
fillOpacity = 0.3,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=5,
fillOpacity = 5)%>%
addMarkers(lng=loc2$long,lat=loc2$lat,
label=~lapply(paste0(loc2[["Location"]],sep="<br>"),HTML),
labelOptions = labelOptions(noHide=T))%>%
addPolygons(data=interim.df%>%
filter(Location=="Greenwich V"),
fillColor = viz.pal[2],
color = "#000000",
weight = 4,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{Location}"))Kebun Baru Mall
leaflet(sg_map%>%
filter(planning_area %in% c("Serangoon","Hougang",
"Sengkang","Ang Mo Kio")))%>%
addProviderTiles("CartoDB.Positron",group = "Greyscale")%>%
addPolygons(fillColor = "white",
color = "grey",
weight = 5,
fillOpacity = 0.3,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=5,
fillOpacity = 5)%>%
addMarkers(lng=loc2$long,lat=loc2$lat,
label=~lapply(paste0(loc2[["Location"]],sep="<br>"),HTML),
labelOptions = labelOptions(noHide=T))%>%
addPolygons(data=interim.df%>%
filter(Location=="Kebun Baru Mall"),
fillColor = viz.pal[3],
color = "#000000",
weight = 4,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{Location}"))We can see that Greenwich Mall has better coverage than the other two locations. So in this blog, we looked at a very easy to follow process of identifying prospective locations for any business.Before we end, I would like to share an important tip that you can use for location icons.
There is a makeIcon function in leaflet package through which you can use the png image of the icons created in power point.Lets see through an example.
Open powerpoint. Then click on insert and then click on icons as shown in the image below.
Once you click on the location, an icon will appear in the slide as shown below.
Now you can color this icon to reflect the color you want, change its size, etc
Once you have the formatting changes, just share it as a png image in your folder
Now you can use this icone shown below
loc_icon<-makeIcon("C:\\Users\\shubh\\OneDrive\\Documents\\R Blogs\\Geospatial\\icon.png",iconWidth = 60,
iconHeight = 60)leaflet(sg_map%>%
filter(planning_area %in% c("Serangoon","Hougang",
"Sengkang","Ang Mo Kio")))%>%
addProviderTiles("CartoDB.Positron",group = "Greyscale")%>%
addPolygons(fillColor = "white",
color = "grey",
weight = 5,
fillOpacity = 0.3,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=5,
fillOpacity = 5)%>%
addMarkers(lng=loc2$long,lat=loc2$lat,
label=~lapply(paste0(loc2[["Location"]],sep="<br>"),HTML),
labelOptions = labelOptions(noHide=T),
icon = loc_icon)%>%
addPolygons(data=interim.df%>%
filter(Location=="Kebun Baru Mall"),
fillColor = viz.pal[3],
color = "#000000",
weight = 4,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{Location}"))We see that icons now have the color we set in the png image.
Step 11: Everything into 1 map
leaflet(sg_map%>%
filter(planning_area %in% c("Serangoon","Hougang",
"Sengkang","Ang Mo Kio")))%>%
addProviderTiles("CartoDB.Positron",group = "Greyscale")%>%
addPolygons(fillColor = "white",
color = "grey",
weight = 5,
fillOpacity = 0.3,
popup = ~glue::glue("<b>{planning_area}</b><br>{district}"))%>%
addPolylines(color="black",
weight=5,
fillOpacity = 5)%>%
addMarkers(lng=loc2$long,lat=loc2$lat,
label=~lapply(paste0(loc2[["Location"]],sep="<br>"),HTML),
labelOptions = labelOptions(noHide=T))%>%
addPolygons(data=interim.df%>%
filter(Location=="NEX Mall Serangoon"),
fillColor = viz.pal[1],
color = "#000000",
weight = 4,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{Location}"),
group = "NEX Mall Serangoon")%>%
addPolygons(data=interim.df%>%
filter(Location=="Greenwich V"),
fillColor = viz.pal[2],
color = "#000000",
weight = 4,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{Location}"),
group = "Greenwich V")%>%
addPolygons(data=interim.df%>%
filter(Location=="Kebun Baru Mall"),
fillColor = viz.pal[3],
color = "#000000",
weight = 4,
fillOpacity = 0.5,
popup = ~glue::glue("<b>{Location}"),
group = "Kebun Baru Mall")%>%
# Layers control allows the user to turn layers on and off
addLayersControl(options = layersControlOptions(collapsed = FALSE),
overlayGroups = c("NEX Mall Serangoon",
"Greenwich V",
"Kebun Baru Mall"))Now you can see that we can use the selection button from the legend pane and visualise the isochrone we want. This helps to compare everything at once.
I hope this blog was useful.In the coming blogs, we will explore more use cases in geospatial.
