# sfdep => for working with sf
# will be using colocation quotients (CLQs), focusing on local (theres also global and pairwise)
::p_load(sf, tmap, tidyverse, sfdep) pacman
In-Class Exercise 5: Advanced Spatial Point Patterns Analysis
Import Data
<- st_read(dsn = "data",
studyArea layer="study_area") %>%
st_transform(crs = 3829) # national projection system of taiwan
Reading layer `study_area' from data source
`/Users/michelle/Desktop/IS415/shelle-mim/IS415-GAA/In-class_Exercise/Wk5/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
<- st_read(dsn = "data",
stores layer="stores") %>%
st_transform(crs = 3829)
Reading layer `stores' from data source
`/Users/michelle/Desktop/IS415/shelle-mim/IS415-GAA/In-class_Exercise/Wk5/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
Visualise the layers
tmap_mode("view")
tm_shape(studyArea) +
tm_polygons()+ #always plot polygon before line
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 Colocation Quotients (LCLQ)
<- include_self(
nb st_knn(st_geometry(stores), 6) # search for the 6 nearest neighbour => each list have 6
# why 6? => so will not have 50-50 (since we include_self, so total 7)
)
<- st_kernel_weights(nb, # calculate weight metrics
wt # target: all stores => convert into a weight metrics
stores, "gaussian",
adaptive = TRUE) # use adaptive method
<- stores %>%
FamilyMart filter(Name == "Family Mart")
<- FamilyMart$Name
A
<- stores %>%
SevenEleven filter(Name == "7-Eleven")
<- SevenEleven$Name
B
# A: target, B: neighbour that we want to find out is colocated or not
<- local_colocation(A, B, nb, wt, 49)
LCLQ
<- cbind(stores, LCLQ) # cannot sort lclq, if not will not match with original data
LCLQ_stores
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))