In-class exercise 5

Published

March 25, 2023

Importing packages

#sfdep - working with simple feature dataframes. For point pattern analysis
#Global colocation quotient 
pacman::p_load(sf, tidyverse, tmap, sfdep)

Importing data

studyArea <- st_read(dsn = 'data',
                     layer ="study_area") %>%
  st_transform(crs = 3829)
Reading layer `study_area' from data source 
  `C:\Study\Y3\S2\IS415\ELAbishek\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
stores <- st_read(dsn = "data",
                  layer = "stores")%>%
  st_transform(crs = 3829)
Reading layer `stores' from data source 
  `C:\Study\Y3\S2\IS415\ELAbishek\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

Visualize the 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))
+−
Name
7-Eleven
Family Mart
Leaflet | Tiles © Esri — Esri, DeLorme, NAVTEQ

Local Colocation Quotients (LCLQ)

nb <- include_self(
  st_knn(st_geometry(stores),6) #6 points, so to avoid a 50-50 situation, we include self so that it is always uneven
)

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

FamilyMart <- stores %>%
  filter(Name == "Family Mart")
A <- FamilyMart$Name

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

#A- target, B- the neighbour we are checking for colocation
LCLQ <- local_colocation(A, B, nb, wt, 49)

LCLQ_stores <- cbind(stores, LCLQ)#Cannot do relational join because the two layers dont ahve unique identifier

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))
+−
X7.Eleven
0.999
Missing
Leaflet | Tiles © Esri — Esri, DeLorme, NAVTEQ
tmap_mode("plot")

Computing kernel weights

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

Preparing the vector list

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

Compute LCLQ

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

Joining output table

LCLQ_stores <- cbind(stores, LCLQ)

Plotting LCLQ values

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))
+−
X7.Eleven
0.999
Missing
Leaflet | Tiles © Esri — Esri, DeLorme, NAVTEQ