In-class Exercise 3: Kernel Density Estimation

Author

Victoria Grace ANN

Published

January 22, 2024

Modified

January 23, 2024

Set up

pacman::p_load(sf, tmap, spNetwork, classInt, virdis, tidyverse, sp)

Importing data

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

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\guacodemoleh\IS415-GAA\In-class_Ex\In-class_Ex03\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

childcare <- st_read(dsn = "data/geospatial", layer="Punggol_CC")

Reading layer `Punggol_CC' from data source 
  `C:\guacodemoleh\IS415-GAA\In-class_Ex\In-class_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM

network <- st_read(dsn = "data/geospatial", layer="Punggol_St")

Reading layer `Punggol_St' from data source 
  `C:\guacodemoleh\IS415-GAA\In-class_Ex\In-class_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM

Examine the structure of the SpatialDataFrame output.

str(network)

Classes 'sf' and 'data.frame':  2642 obs. of  3 variables:
 $ LINK_ID : num  1.16e+08 1.16e+08 1.16e+08 1.16e+08 1.16e+08 ...
 $ ST_NAME : chr  "PUNGGOL RD" "PONGGOL TWENTY-FOURTH AVE" "PONGGOL SEVENTEENTH AVE" "PONGGOL SEVENTEENTH AVE" ...
 $ geometry:sfc_LINESTRING of length 2642; first list element:  'XY' num [1:2, 1:2] 36547 36559 44575 44614
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  ..- attr(*, "names")= chr [1:2] "LINK_ID" "ST_NAME"

str(childcare)

Classes 'sf' and 'data.frame':  61 obs. of  2 variables:
 $ Name    : chr  "kml_10" "kml_99" "kml_100" "kml_101" ...
 $ geometry:sfc_POINT of length 61; first list element:  'XYZ' num  36174 42550 0
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
  ..- attr(*, "names")= chr "Name"

Let’s have a look at how the childcare centres and the network look on the map!

tmap_mode('view') # to zoom
tm_shape(childcare) +
  tm_dots() +
  tm_shape(network) +
  tm_lines()

tmap_mode('plot')

In order to run the Network Constraint Analysis,

Segment the lines
Snap

Creating the line segment,

lixels <- lixelize_lines(network,
                         750, # max dist ppl are willing to walk to a CC
                         mindist=375) #midpoint

samples <- lines_center(lixels)

Performing KDE

More information at https://github.com/JeremyGelb/spNetwork/tree/master

densities <- nkde(network, 
                  events = childcare,
                  w = rep(1,nrow(childcare)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

samples$density <- densities
lixels$density <- densities

# rescaling to help the mapping
samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

Visualise

tmap_mode('view')
tm_shape(lixels)+
  tm_lines(col='density')+
tm_shape(childcare)+
  tm_dots()

tmap_mode('plot')