R Notebook

library(sf)
library(dplyr)
library(spdep)
library(tmap)
library(tidyr)
library(spatstat)

schools  <- st_read("C:/Users/lin92/Downloads/SCHOOL/SCHOOL.shp",options="ENCODING=BIG5",quiet=T)
villages   <- st_read("C:/Users/lin92/Downloads/Taipei_Vill/Taipei_Vill.shp",options="ENCODING=BIG5",quiet=T)
fastfood  <- st_read("C:/Users/lin92/Downloads/Tpe_Fastfood/Tpe_Fastfood.shp",options="ENCODING=BIG5",quiet=T)

PART1-1

fastfood <- st_transform(fastfood, st_crs(villages))
fastfood2 <- st_join(
  fastfood, 
  villages %>% select(TOWN), 
  join = st_within,
  left = FALSE     
)

names(fastfood2)

 [1] "ID"       "CODE"     "ALIAS"    "STORE"    "TYPE_90" 
 [6] "X_COOR"   "Y_COOR"   "TYPE_99"  "TOWN"     "geometry"

unique(fastfood2$TOWNNAME)

NULL

regionA <- c("文山區","大安區","中正區")
regionB <- c("信義區","南港區","松山區")

fastfoodA <- fastfood2 %>% filter(TOWN %in% regionA)
fastfoodB <- fastfood2 %>% filter(TOWN %in% regionB)

nrow(fastfoodA) 

[1] 31

nrow(fastfoodB)

[1] 22

bufA <- st_buffer(fastfoodA, dist = 1000)
bufB <- st_buffer(fastfoodB, dist = 1000)

cntA <- st_intersects(bufA, schools) %>% lengths()
cntB <- st_intersects(bufB, schools) %>% lengths()

ttest_res <- t.test(cntA, cntB, var.equal = FALSE)
print(ttest_res)

    Welch Two Sample t-test

data:  cntA and cntB
t = 1.5977, df = 51, p-value = 0.1163
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1418044  1.2473762
sample estimates:
mean of x mean of y 
 3.870968  3.318182 

P-value>0.05 接受虛無假設，A、B 服務圈內平均學校數無顯著差異

PART1-2

bbox    <- st_bbox(fastfood)
grid500 <- st_make_grid(st_as_sfc(bbox), cellsize = 500, square = TRUE) %>%
           st_sf(grid_id = seq_along(.), geometry = .)

# 計算每格速食店數
ff_cnt  <- st_join(grid500, fastfood, join = st_contains) %>%
           st_drop_geometry() %>%
           count(grid_id) %>%
           rename(count = n)
grid500 <- left_join(grid500, ff_cnt, by = "grid_id") %>%
           mutate(count = replace_na(count, 0))

# Queen 鄰接權重
nb_q    <- poly2nb(grid500, queen = TRUE)
lw_q    <- nb2listw(nb_q, style = "W", zero.policy = TRUE)

# 計算並繪製 Moran’s correlogram
mcorr   <- sp.correlogram(nb_q, var = grid500$count,
                          order = 10, method = "I",
                          style = "W", zero.policy = TRUE)
plot(mcorr, main = "Moran’s Correlogram of Fast-Food Counts")

PART1-3

# 1) 計算 Gi*
Gi_z  <- localG(grid500$count, lw_q)

# 2) Benjamini–Hochberg FDR 校正
pvals <- 2 * pnorm(-abs(Gi_z))
p_adj <- p.adjust(pvals, method = "BH")

# 3) 加回格網並標示顯著熱區
grid500 <- grid500 %>%
  mutate(Gi_z   = as.numeric(Gi_z),
         Gi_p   = pvals,
         Gi_sig = p_adj <= 0.05)

tm_shape(grid500) +
  tm_fill("Gi_z", style = "quantile", title = "Gi* z-score") +
  tm_borders() +
  tm_shape(filter(grid500, Gi_sig)) +
  tm_borders(col = "red", lwd = 2) +
  tm_layout(main.title = "Fast-Food Hot-spots (α=0.05, FDR)")

── tmap v3 code detected ────────────────────────────────────────────
[v3->v4] `tm_fill()`: instead of `style = "quantile"`, use
fill.scale = `tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style' to 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend
of the visual variable `fill` namely 'title' to 'fill.legend =
tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of
`tm_layout(main.title = )`[scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)

PART2-1

library(spatstat.geom) 

win <- as.owin(st_union(villages))
xy <- st_coordinates(schools2)
sch_ppp <- ppp(
  x     = xy[,1],
  y     = xy[,2],
  window= win,
  marks = factor(schools2$TOWN) 
)


levels(marks(sch_ppp))

 [1] "士林區" "大同區" "大安區" "中山區" "中正區" "內湖區" "文山區"
 [8] "北投區" "松山區" "信義區" "南港區" "萬華區"

table(marks(sch_ppp))

士林區 大同區 大安區 中山區 中正區 內湖區 文山區 北投區 松山區 
    20      9     13     11      8     12     22     17      8 
信義區 南港區 萬華區 
     9      7     12 

regions <- c("大安區","文山區","信義區")
F_list  <- lapply(regions, function(reg) {
  Fest(sch_ppp[marks(sch_ppp)==reg], ff_ppp, correction="rs")
})

plot(F_list[[1]], ylim=c(0,1), main="F(d)：各區學校→最近速食店累積分佈", legend=FALSE)
cols <- c("red","blue","darkgreen")
for(i in 2:3) lines(F_list[[i]], col=cols[i])

legend("bottomright", regions, col=cols, lty=1)

NA
NA

PART2-2

library(spatstat.geom)  
library(spatstat) 


r0    <- 3000
K_res <- Kest(ff_ppp, r=seq(0, r0, length=100), correction="Ripley")

idx   <- which.min(abs(K_res$r - r0))

K3000 <- K_res$iso[idx]
L3000 <- sqrt(K3000 / pi)

cat(sprintf("K(3000) = %.3f\nL(3000) = %.3f\n", K3000, L3000))

K(3000) = 67553664.146
L(3000) = 4637.133

env_K <- envelope(
  ff_ppp, 
  fun      = Kest, 
  nsim     = 99, 
  r        = K_res$r,
  correction = "Ripley",
  global   = FALSE,
  savefuns = TRUE
)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68,
69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85,
86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

lo3000 <- env_K$lo[idx]
hi3000 <- env_K$hi[idx]
L_lo   <- sqrt(lo3000 / pi)
L_hi   <- sqrt(hi3000 / pi)

cat(sprintf("K 95%% CI at 3000: [%.3f, %.3f]\n", lo3000, hi3000))

K 95% CI at 3000: [25257089.649, 34347735.988]

cat(sprintf("L 95%% CI at 3000: [%.3f, %.3f]\n", L_lo, L_hi))

L 95% CI at 3000: [2835.416, 3306.543]

plot(env_K, 
     main = "Ripley’s K Envelope\n(fast-food 點圖)",
     legendargs = list(x="topright", legend=c("observed","simulated CI")))
abline(v = r0, col = "grey50", lty = 2)

points(r0, K3000, pch = 16, col = "red")

3

SH_vill_sf <- st_join(schools_sf, villages_sf, join = st_within, left = FALSE)

school_public <- SH_vill_sf %>% filter(Type == "公立")
school_private <- SH_vill_sf %>% filter(Type == "私立")

win <- as.owin(st_bbox(schools_sf))
coords_public <- st_coordinates(st_transform(school_public, 3826))
coords_private <- st_coordinates(st_transform(school_private, 3826))
coords_fastfood <- st_coordinates(st_transform(fastfood_sf, 3826))

pp_public <- ppp(coords_public[,1], coords_public[,2], window = win)
pp_private <- ppp(coords_private[,1], coords_private[,2], window = win)
pp_fastfood <- ppp(coords_fastfood[,1], coords_fastfood[,2], window = win)

警告： data contain duplicated points

F_public <- Fest(pp_public, pp_fastfood, correction = "rs", r = seq(0, 2000, by = 25))
F_private <- Fest(pp_private, pp_fastfood, correction = "rs", r = seq(0, 2000, by = 25))

plot(F_public$r, F_public$rs, type = "l", col = "blue", lwd = 2,
     xlab = "距離 d (公尺)", ylab = "F(d)",
     main = "公立 vs 私立學校到速食店 F(d)", ylim = c(0, 1))  # 加這行！

lines(F_private$r, F_private$rs, col = "red", lwd = 2)

legend("bottomright", legend = c("公立學校", "私立學校"),
       col = c("blue", "red"), lty = 1, lwd = 2)