part3

列出1000 公尺內最高速食店密度的學校

library(sf)
library(spatstat)

# 載入資料
school <- st_read("/Users/harry8995/Downloads/SCHOOL/SCHOOL.shp")

Reading layer `SCHOOL' from data source 
  `/Users/harry8995/Downloads/SCHOOL/SCHOOL.shp' using driver `ESRI Shapefile'
Simple feature collection with 148 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 297078.6 ymin: 2763290 xmax: 312516.7 ymax: 2784542
Projected CRS: TWD97 / TM2 zone 121

fst <- st_read("/Users/harry8995/Downloads/Tpe_Fastfood/Tpe_Fastfood.shp")

Reading layer `Tpe_Fastfood' from data source 
  `/Users/harry8995/Downloads/Tpe_Fastfood/Tpe_Fastfood.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 98 features and 8 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 297198.9 ymin: 2763885 xmax: 312205.7 ymax: 2781148
Projected CRS: TWD97 / TM2 zone 121

tp_vill <- st_read("/Users/harry8995/Downloads/Taipei_Vill/Taipei_Vill.shp")

Reading layer `Taipei_Vill' from data source 
  `/Users/harry8995/Downloads/Taipei_Vill/Taipei_Vill.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 456 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 296094.4 ymin: 2761518 xmax: 317198.9 ymax: 2789180
Projected CRS: TWD97 / TM2 zone 121

# 合併學校與速食店幾何點，取得整體邊界
all_points <- c(st_geometry(school), st_geometry(fst))
all_bounds <- st_bbox(all_points)

# 建立 spatstat 分析視窗
window <- owin(
  xrange = c(all_bounds["xmin"], all_bounds["xmax"]),
  yrange = c(all_bounds["ymin"], all_bounds["ymax"])
)

# 建立學校的 ppp 物件
school_ppp <- ppp(
  x = st_coordinates(school)[,1],
  y = st_coordinates(school)[,2],
  window = window,
  marks = school$Name
)

# 1. Clark-Evans 檢定
# 檢驗學校點分佈是否隨機（假設為完全空間隨機性 CSR）
ce_test <- clarkevans.test(school_ppp, correction = "none", nsim = 999)
print("Clark-Evans Test Result:")

[1] "Clark-Evans Test Result:"

print(ce_test)

    Clark-Evans test
    No edge correction
    Z-test

data:  school_ppp
R = 0.85761, p-value = 0.00092
alternative hypothesis: two-sided

# 2. Ripley’s K 函數
# 計算 K 函數並與隨機分佈比較
k_env <- envelope(school_ppp, Kest, nsim = 99, rank = 1, global = 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.

plot(k_env, main = "Ripley’s K Function for School Distribution", 
     xlab = "Distance (meters)", ylab = "K(r)")

# 原有程式碼：計算速食店數量
fastfood_ppp <- ppp(
  x = st_coordinates(fst)[,1],
  y = st_coordinates(fst)[,2],
  window = window,
  marks = fst$ALIAS
)

# 計算距離矩陣
dist_matrix <- spatstat.geom::crossdist(school_ppp, fastfood_ppp)

# 設定半徑 1000 公尺，計算每間學校速食店數量
radius <- 1000
fastfood_count_1000m <- apply(dist_matrix, 1, function(d) sum(d <= radius))

# 加入學校資料
school$fastfood_count_1000m <- fastfood_count_1000m

# 排序學校依速食店密度（由大到小）
school_sorted <- school[order(-school$fastfood_count_1000m), ]

# 取得前10名的密度門檻值
cutoff_value <- school_sorted$fastfood_count_1000m[10]

# 篩出所有密度 >= cutoff_value 的學校
top_schools <- school_sorted[school_sorted$fastfood_count_1000m >= cutoff_value, ]

# 輸出結果與數量確認
print(as.data.frame(top_schools[, c("Name", "fastfood_count_1000m")]))

NA
NA
NA

輸出top10學校的shapefile

# 縮短欄位名稱以符合 shapefile 限制
names(top_schools)[names(top_schools) == "fastfood_count_1000m"] <- "ff_count"

# 將中文名稱轉換為 ASCII（避免亂碼，選用此方法會替換中文為<xx>格式）
top_schools$Name <- iconv(top_schools$Name, to = "ASCII", sub = "byte")

# 或者，保留中文並指定編碼（視軟體相容性）
# top_schools$Name <- iconv(top_schools$Name, from = "UTF-8", to = "UTF-8")

# 輸出結果與數量確認
print(as.data.frame(top_schools[, c("Name", "ff_count")]))

# 儲存為 shapefile，指定編碼為 UTF-8
st_write(top_schools, "/Users/harry8995/Downloads/SCHOOL/top10_school.shp", append = FALSE, layer_options = "ENCODING=UTF-8")

Deleting layer `top10_school' using driver `ESRI Shapefile'
Writing layer `top10_school' to data source 
  `/Users/harry8995/Downloads/SCHOOL/top10_school.shp' using driver `ESRI Shapefile'
options:        ENCODING=UTF-8 
Writing 16 features with 5 fields and geometry type Point.