Accessibility to COVID-19 vaccination centres in Singapore

Overview

In a bid to combat Covid-19, Singapore has roll out, Starting from January 4 2023, anyone who needs a Covid-19 immunization or a booster shot can come into any vaccination centers for their shots (The Straits Times, 2023). Residents could visit Vaccine Go Where (https://www.gowhere.gov.sg/vaccine) to look for their nearest vaccination center based on their postal code, a web application designed by Government Technology Agency (GovTech) of Singapore. (The Straits Times, 2023)

In this exercise, we attempt to look at the Geospatial properties of accessibility of vaccination centers and how residents might prefer one vaccination center over another using Geospatial analytics.

We load the required library for the exercise below:

…

pacman::p_load(olsrr, corrplot, ggpubr, sf, spdep, GWmodel, tmap, tidyverse, gtsummary, SpatialAcc, ggstatsplot, reshape2)

Geospatial Data Wrangling

Using the URA Master Plan 2014 subzone boundary GIS data, we will load the spatial data with st_read()

mpsz = st_read(dsn="data/geospatial", layer="MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `D:\Allanckw\ISSS624\Hands-on_Ex6\data\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

We then transform the dataset to the SVY21 projection system and ensure that it is valid with the below code

mpsz_svy21 = st_transform(mpsz, 3414)

mpsz_svy21_SubZone = mpsz_svy21

We will then need to create a 250m radius hexagons GIS data. This data set was created by using st_make_grid() of sf package, we will also create an ID for every hexagon.

grid_250  <- st_as_sfc(mpsz_svy21) %>%
  st_make_grid(square = FALSE, cellsize = c(2.5e2, 2.5e2)) %>% 
  st_sf() %>% 
  mutate(id_250 = 1:nrow(.)) # this will be "final" id

We save the hexagon file into a new RDS

saveRDS(grid_250, "data/geospatial/hexagon.rds")

grid_250 = readRDS("data/geospatial/hexagon.rds")
ttm()

tm_shape(grid_250) +
            tmap_options(check.and.fix = TRUE) +
            tm_polygons("id_250", alpha = 0.1) +
            tm_view(set.zoom.limits = c(10,15))

#plot(grid_250)

mpsz_svy21_PLN = tm_shape(mpsz_svy21)+
  tm_polygons("PLN_AREA_N", alpha = 0.1) +

tm_shape(mpsz_svy21_SubZone) +
  tm_polygons("SUBZONE_N", alpha = 0.2)

Public Transport Points

LTA Datamart provides us with the public transport locations, we will attempt to load them here for analysis

busStop = st_read(dsn="data/geospatial/BusStopLocation", layer="BusStop")

Reading layer `BusStop' from data source 
  `D:\Allanckw\ISSS624\Hands-on_Ex6\data\geospatial\BusStopLocation' 
  using driver `ESRI Shapefile'
Simple feature collection with 5160 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

busStop = st_transform(busStop, 3414)

busStop_XY = do.call(rbind, st_geometry(busStop)) %>% 
    as_tibble() %>% setNames(c("X","Y"))

busStop = cbind(busStop, busStop_XY)

The number of Plan Area N is 55, however the default number of categories is only 30, so we will also need to update it by using tmap_options(max.categories = 55) such that all planning zones are shown

Below is the map of the bus stops

tmap_options(max.categories = 55)

mpsz_svy21_tmap = tm_shape(mpsz_svy21)+
  tm_polygons("PLN_AREA_N", alpha = 0.1) +

tm_shape(mpsz_svy21_SubZone) +
  tm_polygons("SUBZONE_N", alpha = 0.2) 
 
mpsz_svy21_tmap +
  tm_shape(busStop) + 
  tm_dots("BUS_STOP_N", legend.show = FALSE) +
  
  tm_view(set.zoom.limits = c(11,15)) 

mrt = st_read(dsn="data/geospatial/TrainStationExit", layer="Train_Station_Exit_Layer")

Reading layer `Train_Station_Exit_Layer' from data source 
  `D:\Allanckw\ISSS624\Hands-on_Ex6\data\geospatial\TrainStationExit' 
  using driver `ESRI Shapefile'
Simple feature collection with 562 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 6134.086 ymin: 27499.7 xmax: 45356.36 ymax: 47865.92
Projected CRS: SVY21

mrt = st_transform(mrt, 3414)

mrt_XY = do.call(rbind, st_geometry(mrt)) %>% 
    as_tibble() %>% setNames(c("X","Y"))

mrt = cbind(mrt, mrt_XY)

Below is the map of the MRT exits

mpsz_svy21_tmap = tm_shape(mpsz_svy21)+
  tm_polygons("PLN_AREA_N", alpha = 0.1) +

tm_shape(mpsz_svy21_SubZone) +
  tm_polygons("SUBZONE_N", alpha = 0.2) 
 
mpsz_svy21_tmap +
  tm_shape(mrt) + 
  tm_dots("stn_name", legend.show = FALSE) +
  
  tm_view(set.zoom.limits = c(11,15)) 

taxi = st_read(dsn="data/geospatial/TaxiStand", layer="TaxiStop")

Reading layer `TaxiStop' from data source 
  `D:\Allanckw\ISSS624\Hands-on_Ex6\data\geospatial\TaxiStand' 
  using driver `ESRI Shapefile'
Simple feature collection with 354 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 6151.732 ymin: 27461.92 xmax: 45468.51 ymax: 48948.45
Projected CRS: SVY21

taxi = st_transform(taxi, 3414)

taxi_XY = do.call(rbind, st_geometry(taxi)) %>% 
    as_tibble() %>% setNames(c("X","Y"))

taxi = cbind(taxi, taxi_XY)

Below is the map of the taxi stands

mpsz_svy21_tmap = tm_shape(mpsz_svy21)+
  tm_polygons("PLN_AREA_N", alpha = 0.1) +

tm_shape(mpsz_svy21_SubZone) +
  tm_polygons("SUBZONE_N", alpha = 0.2) 
 
mpsz_svy21_tmap +
  tm_shape(taxi) + 
  tm_dots(legend.show = FALSE) +
  
  tm_view(set.zoom.limits = c(11,15)) 

Aspatial Data Wrangling

vacc_centers = read_csv("data/aspatial/vacc_center_cleaned.csv") 

We can use st_as_sfto create a dataframe from the longitude (x) and latitude (y) values. The EPSG 4326 code is used as the dataset is referencing WGS84 geographic coordinate system. We could use st_crs()to verify the coordinate system from the object.

vacc_centers_sf = st_as_sf(vacc_centers, coords = c("Longitude", "Latitude"), crs=4326)
st_crs(vacc_centers_sf)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

vacc_centers_sf = st_transform(vacc_centers_sf, 3414)
saveRDS(vacc_centers_sf, "data/aspatial/vacc_centers_sf.rds")

vacc_centers_sf = readRDS("data/aspatial/vacc_centers_sf.rds")

Below is the map of the vaccination centers

mpsz_svy21_tmap = tm_shape(mpsz_svy21)+
  tm_polygons("PLN_AREA_N", alpha = 0.1) +

tm_shape(mpsz_svy21_SubZone) +
  tm_polygons("SUBZONE_N", alpha = 0.2) 
 
mpsz_svy21_tmap +
  tm_shape(vacc_centers_sf) + 
  tm_dots("name", legend.show = FALSE) +
  
  tm_view(set.zoom.limits = c(11,15)) 

Reference

The Straits Times (2023), All can walk in for Covid-19 jabs at vaccination centres from Jan 4 https://www.straitstimes.com/singapore/people-can-walk-into-vaccine-centres-to-get-covid-19-jab-from-wednesday