In-class Ex 5: Advanced Spatial Point Pattern Analysis

Published

February 6, 2023

Modified

February 12, 2023

Installing packages

sfdep: spatial and spatial time cubes, colocation purposes

pacman::p_load(tidyverse, sf, tmap, sfdep)

Importing Data

Transform to Taiwan’s projection System (3829)

studyArea <- st_read(dsn = "data",
                    layer = "study_area") %>%
  st_transform(crs = 3829)

Reading layer `study_area' from data source 
  `C:\zoe-chia\IS415\In-class_Ex\In-class_Ex05\data' using driver `ESRI Shapefile'
Simple feature collection with 7 features and 7 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 121.4836 ymin: 25.00776 xmax: 121.592 ymax: 25.09288
Geodetic CRS:  TWD97

stores <- st_read(dsn = "data", 
                  layer = "stores") %>%
  st_transform(crs = 3829)

Reading layer `stores' from data source 
  `C:\zoe-chia\IS415\In-class_Ex\In-class_Ex05\data' using driver `ESRI Shapefile'
Simple feature collection with 1409 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 121.4902 ymin: 25.01257 xmax: 121.5874 ymax: 25.08557
Geodetic CRS:  TWD97

Visualising the sf layers

tmap_mode('view')
tm_shape(studyArea) + 
  tm_polygons() + 
  tm_shape(stores) + 
  tm_dots(col = "Name", 
          size = 0.01, 
          border.col = "black", 
          border.lwd = 0.5) + 
  tm_view(set.zoom.limits = c(12, 16))

Local Coloation Quotients (LCLQ)

Calculating nearest neighbour list:

nb <- include_self(
  st_knn(st_geometry(stores), 6))
# 6 means to search for the 6 nearest neighbours
# - Adaptive method
# - Total 7 points (including self) so that you can get a good split. Avoid even splits e.g. 3, 3

wt <- st_kernel_weights(nb, 
                        stores, 
                        "gaussian",
                        adaptive = TRUE
                        )

FamilyMart <- stores %>%
  filter(Name == "Family Mart") # Select stores where name = 'Family Mart', in dataframe form 
A <- FamilyMart$Name #vector, target

SevenEleven <- stores %>%
  filter(Name == "7-Eleven")
B <- SevenEleven$Name

Run 50 simulations:

LCLQ <- local_colocation(A, B, nb, wt, 49)

Combine stores and local location table. - Not doing left join / right join because these two layers have no unique identifiers. - Don’t sort data by code - Left hand side should be the geomerty, ‘stores’

LCLQ_stores <- cbind(stores, LCLQ)

tmap_mode('view')
tm_shape(studyArea) + 
tm_polygons() + 
  tm_shape(LCLQ_stores) + 
  tm_dots(col = "X7.Eleven", 
          size = 0.01, 
          border.col = "black", 
          border.lwd = 0.5) + 
  tm_view(set.zoom.limits = c(12, 16))

tmap_mode('plot')