Tuesday, July 30, 2024

Geospatial Tutorial 2 - Drive Time Analysis in Singapore

Geospatial Analysis Tutorial 2: Identifying locations for Opening Restaurents

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.3696


loc2<-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     3


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))

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.df
Simple 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.

at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)

Price Elasticity Model in Python

Price Elasticity Analysis in Python Introduction ¶ In this blog we will look at how to find the price e...