In-class Exercise 4

In this in class exercise, I gain a deeper understanding on the different functions to perform spatial point patterns analysis.

1. Installing and loading required packages

My notes:

• sf is mainly used for importing geospatial data
• raster is for handling

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)
}
# onemapsgapi

1.1 More information about the packages used

My notes:

• maptools: Manipulate geographic data.

  • Will be used to convert spatial objects into ppp format of spatstat

• sf: used for importing, managing, and processing geospatial data

  • specifically vector-based geospatial data

• raster: reads, writes, manipulates, analyses and model of gridded spatial data

• spatstat: contains useful functions for point pattern analysis.

  • Will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.

• tmap: provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.

• tidyverse (for CHILDCARE dataset): used for importing, wrangling and visualising data and consists of a family of packages:

  • readr for importing csv data,
  • readxl for importing Excel worksheet,
  • tidyr for manipulating data,
  • dplyr for transforming data, and
  • ggplot2 for visualising data

2. Import Geospatial Data

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

Reading layer `CostalOutline' from data source 
  `C:\aisyahajit2018\IS415\IS415_blog\_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:\aisyahajit2018\IS415\IS415_blog\_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

3. Importing Aspatial Data

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

My notes:

** Difference between readRDS and read_rds:

• discouraged from using readRDS (baseR method) to import data as they are relatively slow when dealing with big data
• they are not tibble data format
• Hence, better to use read_rds

3.1 Converting from aspatial to geospatial

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

My notes:

• For this apatial data, we have 2 coordinate data one is in degree and one is in metres
• as we know that the data is by SLA, so we used

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)

4. Geospatial Data Wrangling

• st_set_crs is for assigning crs with the EPSG code
• Difference between SVY21 and WGS 84:
  • Allow us to know the data itself
  • Originally, the CHILDCARE is in latlng and is measured in national degree.
  • EPSG code for WGS8 is 4326
  • 1st example, we already know it is in metres, therefore we straightaway assign the crs to 3414.
  • For 2nd example, since it is latlng, we convert it to 4326 (wgs 84 format) first then transform the crs.
• Difference between set and transform:
  • When we set using st_set_crs, we only assign and not performing any mathematical transformation.
  • But when we use st_transform, we are actually converting decimal degree to metres.
• NOTE: From sf to spatial* Classes to sp format, there will some sort of information lost

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

4.2 Convert from Spatial* to sp format

• Will have the projection and coordinates

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

My notes:

• Difference between spatial df vs spatial class:
  • When working with sp layer, there is no associated table. So all the attribute information will be gone.
  • We will only be left with the coordinates etc

4.3 Convert from sp to spatstat ppp format

• Provide you with the explicit data model (e.g owin/rectangle)

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

plot(childcare_ppp)

plot(CHAS_ppp)

summary(childcare_ppp)

Planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units

summary(CHAS_ppp)

Planar point pattern:  1192 points
Average intensity 5.39042e-07 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: rectangle = [0, 45475.65] x [0, 48626.7] units
Window area = 2211330000 square units

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)

Here we are switching back to plot mode:

tmap_mode('plot')