In-Class Exercise 4

Today’s in-class exercise will bring us through a step-by-step of data preparation + wrangling, as well as letting us test out the different functions we’re learned for analysis so far ❤️

1.0 Overview

In Hands-On Exercise 4, we learned how to perform Spatial Point Patterns Analysis with spatstat, SPPA for short! As such, in this In-Class Exercise, we’ll be putting what we learned into use!

2.0 Setup

2.1 Packages Used

The R packages we’ll be introducing today are:

• spatstat: used for point pattern analysis
• raster: reads, writes, manipulates, analyses and models gridded spatial data (i.e. raster-based geographical data)
• maptools: a set of tools for manipulating geographic data

In addition, we’ll be using the packages from our last exercise:

• sf: used for importing, managing, and processing geospatial data
• tmap: used for creating thematic maps, such as choropleth and bubble maps
• tidyverse: used for importing, wrangling and visualising data (and other data science tasks!)

packages = c('maptools', 'sf', 'raster','spatstat', 'tmap', 'tidyverse')
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p,character.only = T)
}

2.2 Data Used

The datasets used for this exercise are:

• CHILDCARE, a point feature data in rds format from data.gov.sg, providing both location and attribute information of childcare centres
• CHAS, a point feature data in rds format that provides location and attribute information of CHAS clinics
• MP14_SUBZONE_WEB_PL, a polygon feature data in ESRI shapefile format from data.gov.sg, providing information of URA 2014 Master Plan Planning Subzone boundary data
• CostalOutline, a polygon feature data in ESRI shapefile format from sla.gov.sg, showing the national boundary of Singapore

2.3 Importing Data

2.3.1 Importing Geospatial Data

sg_sf <- st_read(dsn = "data/shapefile", 
                 layer="CostalOutline")

Reading layer `CostalOutline' from data source 
  `C:\IS415\IS415_msty\_posts\2021-09-06-in-class-exercise-4\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21

mpsz_sf <- st_read(dsn = "data/shapefile", 
                layer = "MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\IS415\IS415_msty\_posts\2021-09-06-in-class-exercise-4\data\shapefile' 
  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

2.3.2 Importing Aspatial Data

childcare <- read_rds("data/rds/childcare.rds")
CHAS <- read_rds("data/rds/CHAS.rds")

2.4 Data Preparation + Wrangling

2.4.1 Converting from Aspatial to Geospatial

Callback to Hands-On Exercise 02 - converting from aspatial data to geospatial data with our st_as_sf() function 💪

CHAS_sf <- st_as_sf(CHAS,
                    coords = c("X_COORDINATE",
                               "Y_COORDINATE"),
                    crs=3414)

childcare$Lat <- as.numeric(childcare$Lat)
childcare$Lng <- as.numeric(childcare$Lng)
childcare_sf <- st_as_sf(childcare,
                    coords = c("Lng",
                               "Lat"),
                    crs=4326) %>%
  st_transform(crs=3414)

Callback: why do we set crs at 3414 at the end? When converting aspatial data into geospatial, we usually used particular ‘geograhpical’ variables (x-coord and y-coord, Lat and Lng) to situate the data - but we also need to know the coordinate system being used. Since our data is from government sources, and we know that the SLA tends to use the Singapore Projected Coordinate system (based on official documentation), We should have our shapefiles in SVY21 format - whose EPSG code is 3414.

💡 Thinking Point: why do we immediately crs as 3414 for CHAS, but assign as 4326 and then transform it to 3414 for childcare? Because childcare’s ‘geographical’ variables are longtitude and latitde, in demical degrees. The most appropriate coordinate system for that is WGS84, crs code 4326 - and then transform it to our desired crs 3414 👍

2.4.2 Converting from sf to Spatial* classes

childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)

2.4.3 Converting from Spatial* classes to sp format

childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS,"SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")

2.4.4 Converting from sp format to spatstat ppp format

childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")

3.0 Analysis

Our goal is to compare the distribution of childcare centre vs CHAS clinics. They both have relatively large data points and are rather scattered, making them ideal to compare, and it’d be interesting to see if there are any correlations between them as their target groups are largely different: children vs senior citizens.

3.1 Interactive Overview

tmap_mode('view')
tm_shape(childcare_sf) +
  tm_dots(alpha=0.4,
          col="blue",
          size=0.05) +
tm_shape(CHAS_sf) +
  tm_dots(alpha=0.4,
          col="red",
          size=0.05)