空間分析 Final Report

#讀取資料

health_sf = st_read("./data/臺北市健康餐飲.kml", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `臺北市健康餐飲' from data source 
  `D:\R\spatial analysis\Final\data\臺北市健康餐飲.kml' using driver `KML'
Simple feature collection with 234 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: 121.4936 ymin: 24.96722 xmax: 121.6181 ymax: 25.13751
Geodetic CRS:  WGS 84

health_sf <- st_transform(health_sf, crs = 3826)
SCHOOL_sf = st_read("./data/SCHOOL.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `SCHOOL' from data source `D:\R\spatial analysis\Final\data\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

Taipei_Vill_sf = st_read("./data/Taipei_Vill.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Taipei_Vill' from data source `D:\R\spatial analysis\Final\data\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

Tpe_Fastfood_sf = st_read("./data/Tpe_Fastfood.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Tpe_Fastfood' from data source `D:\R\spatial analysis\Final\data\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

Part 1 - 實作題

第一題

假設速食店的服務範圍是 1 公里，比較 A 區(文山+大安+中正) 與 B 區(信義+南港+松山) 這兩個地區的每一家速食店在服務可及範圍內，所涵蓋學校數量的平均值，是否有統計顯著差異。(需列出虛無假設與對立假設，統計檢定量，以及檢定的顯著水準等)。

• \(H_0\): A區平均值等於B區平均值
• \(H_A\): A區平均值不等於B區平均值
• 顯著水準：\(\alpha = 0.05\)

FAST_COUNTY_sf = st_join(Tpe_Fastfood_sf, Taipei_Vill_sf, join = st_within)
FAST_buffer_sf = st_buffer(FAST_COUNTY_sf, dist = 1000)
FAST_SCHOOL_sf = st_join(SCHOOL_sf, FAST_buffer_sf, join = st_within)

FAST_A_sf = filter(FAST_SCHOOL_sf, TOWN %in% c("文山區", "大安區", "中正區"))
FAST_B_sf = filter(FAST_SCHOOL_sf, TOWN %in% c("信義區", "南港區", "松山區"))

count_A_table = table(FAST_A_sf$ID)
count_B_table = table(FAST_B_sf$ID)

count_A_df = as.data.frame(count_A_table)
colnames(count_A_df) <- c("ALIAS","學校數")
count_B_df = as.data.frame(count_B_table)
colnames(count_B_df) <- c("ALIAS","學校數")

t.test(count_A_df$學校數, count_B_df$學校數, alternative = "two.sided")

    Welch Two Sample t-test

data:  count_A_df$學校數 and count_B_df$學校數
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.1163>\(\alpha\)，所以無法拒絕虛無假設，因此兩區的涵蓋學校數量的平均值，沒有統計顯著差異

第二題

# 建立 500 公尺網格
grid500 <- st_make_grid(Taipei_Vill_sf, cellsize = 500, square = TRUE) %>% 
  st_sf() %>% 
  st_intersection(st_union(Taipei_Vill_sf)) %>% 
  st_sf()
st_crs(grid500) <- 3826

# 計算每格內速食店數量
grid500$Fastfood_N <- lengths(st_intersects(grid500, Tpe_Fastfood_sf))

# 建立 contiguity 鄰接關係 (queen)
nb_queen <- poly2nb(grid500, queen = TRUE)
lw_queen <- nb2listw(nb_queen, style = "W")

# 繪製網格 + 速食店數量
tmap_mode("plot")
tm_shape(grid500) +
  tm_polygons("Fastfood_N", palette = "YlOrRd", title = "the amount of fast food restaurants") +
  tm_layout(main.title = "500-meter grid of fast food restaurants")

# 計算 Moran’s I Correlogram（order = 10）
moran_corr <- sp.correlogram(nb_queen, grid500$Fastfood_N, order = 10, method = "I", style = "W")

# 繪圖
plot(moran_corr, main = "Moran's I Correlogram of Fast Food Counts")

##### 額外：使用FDR校正 #####

# 計算前 10 階的 Moran's I
moran_corr <- sp.correlogram(nb_q, grid500$Fastfood_N, order = 10, method = "I", style = "W", zero.policy = TRUE)

# 轉成 dataframe 繪圖用
df_corr <- data.frame(
  Order = 1:10,
  I = moran_corr$res[, 1],  # 第一欄是 I
  p = moran_corr$res[, 2]   # 第二欄是 p 值
)

# 繪圖
ggplot(df_corr, aes(x = Order, y = I)) +
  geom_line() +
  geom_point(aes(color = p < 0.05), size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
  scale_color_manual(values = c("FALSE" = "grey", "TRUE" = "red"), labels = c("Not Significant", "Significant")) +
  labs(title = "Moran's I Correlogram of Fast Food Counts",
       subtitle = "Spatial autocorrelation by contiguity order",
       x = "Contiguity Order", y = "Moran's I", color = "p < 0.05") +
  theme_minimal()

##### 額外：使用FDR校正 #####

# 建立 500m 網格
grid500 <- st_make_grid(Taipei_Vill_sf, cellsize = 500, square = TRUE) %>% 
  st_sf() %>%
  st_intersection(st_union(Taipei_Vill_sf)) %>%
  st_sf()
st_crs(grid500) <- 3826

# 加入每格速食店數量
grid500$Fastfood_N <- lengths(st_intersects(grid500, Tpe_Fastfood_sf))

# 建立鄰接關係（Queen）
nb_q <- poly2nb(grid500, queen = TRUE)
lw_q <- nb2listw(nb_q, style = "W", zero.policy = TRUE)

# 計算 Local Moran's I
lisa <- localmoran(grid500$Fastfood_N, lw_q, zero.policy = TRUE)

# 加入結果並 FDR 校正
grid500$Ii         <- lisa[,1]
grid500$EIi        <- lisa[,2]
grid500$VarIi      <- lisa[,3]
grid500$Z.Ii       <- lisa[,4]
grid500$p.Ii       <- lisa[,5]
grid500$p.Ii.fdr   <- p.adjust(lisa[,5], method = "fdr")

# LISA 分類（依照顯著性 + 正負相關）
lag_FF <- lag.listw(lw_q, grid500$Fastfood_N)
grid500$lag_FF <- lag_FF

grid500$LISA_type <- "Not significant"
grid500$LISA_type[grid500$Fastfood_N > mean(grid500$Fastfood_N) & 
                  grid500$lag_FF > mean(grid500$lag_FF) & 
                  grid500$p.Ii.fdr < 0.05] <- "High-High"
grid500$LISA_type[grid500$Fastfood_N < mean(grid500$Fastfood_N) & 
                  grid500$lag_FF < mean(grid500$lag_FF) & 
                  grid500$p.Ii.fdr < 0.05] <- "Low-Low"
grid500$LISA_type[grid500$Fastfood_N > mean(grid500$Fastfood_N) & 
                  grid500$lag_FF < mean(grid500$lag_FF) & 
                  grid500$p.Ii.fdr < 0.05] <- "High-Low"
grid500$LISA_type[grid500$Fastfood_N < mean(grid500$Fastfood_N) & 
                  grid500$lag_FF > mean(grid500$lag_FF) & 
                  grid500$p.Ii.fdr < 0.05] <- "Low-High"

grid500$LISA_type <- factor(grid500$LISA_type,
                            levels = c("High-High", "Low-Low", "High-Low", "Low-High", "Not significant"))

# 畫出 LISA 類型地圖
tmap_mode("plot")

ℹ tmap mode set to "plot".

tm_shape(grid500) +
  tm_polygons("LISA_type", palette = c("red", "blue", "orange", "skyblue", "grey80"),
              title = "LISA Type\n(FDR Corrected)") +
  tm_layout(main.title = "Local Moran's I for Fast Food Clusters in Taipei",
            legend.outside = TRUE)

── tmap v3 code detected ──────────────────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of the visual variable `fill` namely
'palette' (rename to 'values') to fill.scale = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: 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 = )`

速食店在空間上具有顯著的聚集現象，尤其在第 1 階（500 公尺鄰近）時 Moran’s I 約為 0.18，顯著高於 0，代表空間自相關強烈。

隨著鄰接階數增加，Moran’s I 值逐漸下降，但在前 6 階（即約 3 公里範圍內）仍維持正向且有顯著性，表示聚集效應具有一定的空間延伸範圍。

自第 7 階之後，Moran’s I 值接近 0，且信賴區間與 0 重疊，顯示自相關趨於隨機，空間聚集效應已不明顯。

總結： 速食店在 500 公尺至 3 公里範圍內存在顯著空間聚集現象，顯示店家選址或商圈形成具地理集中性；超過此距離後，自相關趨近隨機，表示空間影響效應已消散。

第三題

# --- STEP 0: 建立學校數量欄位（與前面一樣） ---
grid500$School_N <- lengths(st_intersects(grid500, SCHOOL_sf))

# 建立 Queen 鄰接權重（如果還沒執行過）
nb_q <- poly2nb(grid500, queen = TRUE)
lw_q <- nb2listw(nb_q, style = "W", zero.policy = TRUE)

# 找出有鄰居的格子
valid_idx <- which(card(nb_q) > 0)

# 初始化並計算 Gi* 統計量
Gi_Z_FF <- rep(NA, nrow(grid500))
Gi_Z_FF[valid_idx] <- localG(grid500$Fastfood_N[valid_idx], lw_q, zero.policy = TRUE)

grid500$Gi_Z_FF <- Gi_Z_FF
grid500$p_FF <- 2 * pnorm(-abs(grid500$Gi_Z_FF))  # 雙尾 p 值
grid500$p_FF_fdr <- p.adjust(grid500$p_FF, method = "fdr")  # FDR 修正
grid500$Hotspot_FF <- ifelse(grid500$p_FF_fdr < 0.05 & grid500$Gi_Z_FF > 0, "Hotspot", "Not significant")

# 1. 計算 Gi* Z 值（學校）
Gi_Z_SC <- rep(NA, nrow(grid500))
Gi_Z_SC[valid_idx] <- localG(grid500$School_N[valid_idx], lw_q, zero.policy = TRUE)
grid500$Gi_Z_SC <- Gi_Z_SC

# 2. 計算 p 值與 FDR 校正
grid500$p_SC <- 2 * pnorm(-abs(grid500$Gi_Z_SC))
grid500$p_SC_fdr <- p.adjust(grid500$p_SC, method = "fdr")

# 3. 正確產出 Hotspot_SC
grid500$Hotspot_SC <- ifelse(
  is.na(grid500$p_SC_fdr), "Not significant",
  ifelse(grid500$p_SC_fdr < 0.05 & grid500$Gi_Z_SC > 0, "Hotspot", "Not significant")
)

# 4. 設定為 factor，避免 tmap 判定為 NA
grid500$Hotspot_SC <- factor(grid500$Hotspot_SC, levels = c("Hotspot", "Not significant"))
# --- STEP 2: 熱區圖繪製 ---

# 速食店熱區地圖
tm_shape(grid500) +
  tm_polygons("Hotspot_FF", palette = c("red", "grey80"),
              title = "Fast Food Gi* Hotspot") +
  tm_layout(main.title = "Fast Food Hotspots (Gi*, FDR corrected)",
            legend.outside = TRUE)

# 學校熱區地圖
tm_shape(grid500) +
  tm_polygons("Hotspot_SC", palette = c("blue", "grey80"),
              title = "School Gi* Hotspot") +
  tm_layout(main.title = "School Hotspots (Gi*, FDR corrected)",
            legend.outside = TRUE)

結論（速食店與學校熱區）： 速食店熱區（紅色）主要集中於圖中南側偏中與東南區域，呈現連續聚集，顯示這些區域的速食店密度顯著高於其他地區。 學校熱區（藍色）則較為零星分布，雖也多集中於南側，但與速食店的熱區僅部分重疊，整體空間分布差異明顯。 兩者空間分布顯示：速食店熱區與學校熱區不完全一致，但有重疊。

Part 2

第一題

#讀取資料

fastfood <- st_read("./data/Tpe_Fastfood.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Tpe_Fastfood' from data source `D:\R\spatial analysis\Final\data\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

villages <- st_read("./data/Taipei_Vill.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Taipei_Vill' from data source `D:\R\spatial analysis\Final\data\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

schools <- st_read("./data/SCHOOL.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `SCHOOL' from data source `D:\R\spatial analysis\Final\data\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

#Step 2：選出三個行政區的學校（大安、文山、信義）
villages <- villages %>%
  filter(grepl("大安區|文山區|信義區", TOWN))  # 注意欄位名稱可能需對應你的資料，如 DISTRICT_NAME

#空間相交找出屬於三區的學校
schools_three <- st_join(schools, villages, join = st_within) %>%
  filter(!is.na(TOWN))  # 確保有交集成功

#計算每所學校到所有速食店的距離，並做 F(d) 函數

#定義距離門檻（以公尺為單位），這裡設定從 0 到 1000 公尺，每 100 公尺一個門檻
d_values <- seq(0, 1000, by = 100)

#建立一個函數計算 F(d) 值
calc_Fd <- function(schools_sf, fastfood_sf, d_seq) {
  result <- data.frame(distance = d_seq, Fd = NA)
  n_school <- nrow(schools_sf)
  
  for (i in seq_along(d_seq)) {
    d <- set_units(d_seq[i], "m") #單位為公尺
    
    #每間學校是否在 d 公尺內至少有一間速食店？
    count_within_d <- sum(sapply(1:n_school, function(j) {
      any(st_distance(schools_sf[j,], fastfood_sf) <= d)
    }))
    
    #F(d) = 滿足條件的學校數 / 總學校數
    result$Fd[i] <- count_within_d / n_school
  }
  return(result)
}

#分三區做 F(d)
fd_da_an <- calc_Fd(schools_three %>% filter(grepl("大安區", TOWN)), fastfood, d_values)
fd_wen_shan <- calc_Fd(schools_three %>% filter(grepl("文山區", TOWN)), fastfood, d_values)
fd_xin_yi <- calc_Fd(schools_three %>% filter(grepl("信義區", TOWN)), fastfood, d_values)

#加上區名
fd_da_an$District <- "大安區"
fd_wen_shan$District <- "文山區"
fd_xin_yi$District <- "信義區"

#合併資料
fd_all <- rbind(fd_da_an, fd_wen_shan, fd_xin_yi)

#繪圖比較 F(d)
ggplot(fd_all, aes(x = distance, y = Fd, color = District)) +
  geom_line(size = 1.2) +
  scale_x_continuous(breaks = d_values) +
  labs(
    title = "F(d) 函數比較：三區學校周邊速食店分布",
    x = "距離 d（公尺）",
    y = "F(d) = 至少一間速食店的學校比例"
  ) +
  theme_minimal() +
  theme(legend.title = element_blank())

NA
NA

set.seed(42)  #為了結果可重現

#只處理大安區
schools_da_an <- schools_three %>% filter(grepl("大安區", TOWN))
n_schools <- nrow(schools_da_an)
da_an_area <- villages %>% filter(TOWN == "大安區")

#模擬次數
n_sim <- 99

#建立一個空的 data frame 儲存模擬結果
sim_Fd_all <- matrix(NA, nrow = length(d_values), ncol = n_sim)

#進行模擬
for (i in 1:n_sim) {
  # 在大安區隨機產生 n_schools 個點
  random_pts <- st_sample(da_an_area, size = n_schools, type = "random") %>%
    st_as_sf()
  
  # 計算 F(d)
  sim_fd <- calc_Fd(random_pts, fastfood, d_values)
  
  # 存入模擬結果矩陣
  sim_Fd_all[, i] <- sim_fd$Fd
}

#產生模擬的 envelopes（最大值與最小值）
Fd_sim_min <- apply(sim_Fd_all, 1, min)
Fd_sim_max <- apply(sim_Fd_all, 1, max)

#合併為一個資料框
fd_envelope <- data.frame(
  distance = d_values,
  min = Fd_sim_min,
  max = Fd_sim_max
)

#加上實際的 F(d)
fd_actual <- fd_da_an  # 剛剛計算過的大安區實際 F(d)

#Step 6：繪圖（實際 vs. 模擬區間）
ggplot() +
  geom_ribbon(data = fd_envelope, aes(x = distance, ymin = min, ymax = max), 
              fill = "lightblue", alpha = 0.4) +
  geom_line(data = fd_actual, aes(x = distance, y = Fd), color = "red", size = 1.2) +
  labs(
    title = "大安區學校周邊速食店的 F(d) 與模擬區間",
    subtitle = "紅線為實際值，藍色區域為模擬 99 次所得之區間",
    x = "距離 d（公尺）",
    y = "F(d)"
  ) +
  theme_minimal()

set.seed(42)

#只處理信義區
schools_da_an <- schools_three %>% filter(grepl("信義區", TOWN))
n_schools <- nrow(schools_da_an)
da_an_area <- villages %>% filter(TOWN == "信義區")

#模擬次數
n_sim <- 99

#建立一個空的 data frame 儲存模擬結果
sim_Fd_all <- matrix(NA, nrow = length(d_values), ncol = n_sim)

#進行模擬
for (i in 1:n_sim) {
  #信義區隨機生成 n_schools 個點
  random_pts <- st_sample(da_an_area, size = n_schools, type = "random") %>%
    st_as_sf()
  
  # 計算 F(d)
  sim_fd <- calc_Fd(random_pts, fastfood, d_values)
  
  # 存入模擬結果矩陣
  sim_Fd_all[, i] <- sim_fd$Fd
}

#產生模擬的 envelopes（最大值與最小值）
Fd_sim_min <- apply(sim_Fd_all, 1, min)
Fd_sim_max <- apply(sim_Fd_all, 1, max)

#合併為一個資料框
fd_envelope <- data.frame(
  distance = d_values,
  min = Fd_sim_min,
  max = Fd_sim_max
)

#加上實際的 F(d)
fd_actual <- fd_da_an  # 剛剛計算過的信義區實際 F(d)

#Step 6：繪圖（實際 vs. 模擬區間）
ggplot() +
  geom_ribbon(data = fd_envelope, aes(x = distance, ymin = min, ymax = max), 
              fill = "lightblue", alpha = 0.4) +
  geom_line(data = fd_actual, aes(x = distance, y = Fd), color = "red", size = 1.2) +
  labs(
    title = "信義區學校周邊速食店的 F(d) 與模擬區間",
    subtitle = "紅線為實際值，藍色區域為模擬 99 次所得之區間",
    x = "距離 d（公尺）",
    y = "F(d)"
  ) +
  theme_minimal()

set.seed(42)

#只處理文山區
schools_da_an <- schools_three %>% filter(grepl("文山區", TOWN))
n_schools <- nrow(schools_da_an)
da_an_area <- villages %>% filter(TOWN == "文山區")

#模擬次數
n_sim <- 99

#建立一個空的 data frame 儲存模擬結果
sim_Fd_all <- matrix(NA, nrow = length(d_values), ncol = n_sim)

#進行模擬
for (i in 1:n_sim) {
  #文山區隨機生成 n_schools 個點
  random_pts <- st_sample(da_an_area, size = n_schools, type = "random") %>%
    st_as_sf()
  
  # 計算 F(d)
  sim_fd <- calc_Fd(random_pts, fastfood, d_values)
  
  # 存入模擬結果矩陣
  sim_Fd_all[, i] <- sim_fd$Fd
}

#產生模擬的 envelopes（最大值與最小值）
Fd_sim_min <- apply(sim_Fd_all, 1, min)
Fd_sim_max <- apply(sim_Fd_all, 1, max)

#合併為一個資料框
fd_envelope <- data.frame(
  distance = d_values,
  min = Fd_sim_min,
  max = Fd_sim_max
)

#加上實際的 F(d)
fd_actual <- fd_da_an  # 剛剛計算過的文山區實際 F(d)

#Step 6：繪圖（實際 vs. 模擬區間）
ggplot() +
  geom_ribbon(data = fd_envelope, aes(x = distance, ymin = min, ymax = max), 
              fill = "lightblue", alpha = 0.4) +
  geom_line(data = fd_actual, aes(x = distance, y = Fd), color = "red", size = 1.2) +
  labs(
    title = "文山區學校周邊速食店的 F(d) 與模擬區間",
    subtitle = "紅線為實際值，藍色區域為模擬 99 次所得之區間",
    x = "距離 d（公尺）",
    y = "F(d)"
  ) +
  theme_minimal()

第二題

win_polygon <- st_as_sfc(st_bbox(Tpe_Fastfood_sf))
win <- as.owin(win_polygon)

coords <- st_coordinates(Tpe_Fastfood_sf)
fastfood_ppp <- ppp(x = coords[, 1], y = coords[, 2], window = win)

r_seq <- seq(0, 5000, by = 100)

k_result <- Kest(fastfood_ppp, r = r_seq)
k_3000 <- k_result$iso[which(k_result$r == 3000)]
l_3000 <- sqrt(k_3000 / pi) -3000

cat("K(3000) =", k_3000, "\n")

K(3000) = 60111493 

cat("L(3000) =", l_3000, "\n")

L(3000) = 1374.252 

env <- envelope(fastfood_ppp,
                fun = Kest,
                nsim = 99,
                r = r_seq,
                correction = "iso",
                savefuns = TRUE,
                savepatterns = TRUE,
                verbose = FALSE)

r_index <- which(env$r == 3000)
upper_L <- sqrt(env$hi[r_index] / pi) -3000
lower_L <- sqrt(env$lo[r_index] / pi) -3000

cat("95% CI for L(3000): [", lower_L, ",", upper_L, "]\n")

95% CI for L(3000): [ -159.3052 , 162.5662 ]

plot(env, . ~ r, main = "Ripley's k-function with 95% CI")
abline(v = 3000, col = "red", lty = 2)

第三題

public_sf <- subset(SCHOOL_sf, Type == "公立")
private_sf <- subset(SCHOOL_sf, Type == "私立")

coords_pub <- st_coordinates(public_sf)
coords_pri <- st_coordinates(private_sf)
coords_ff  <- st_coordinates(Tpe_Fastfood_sf)

win <- as.owin(st_as_sfc(st_bbox(SCHOOL_sf)))

pp_all_pub <- ppp(x = c(coords_pub[, 1], coords_ff[, 1]),
                  y = c(coords_pub[, 2], coords_ff[, 2]),
                  window = win,
                  marks = factor(c(rep("school", nrow(coords_pub)),
                                   rep("fastfood", nrow(coords_ff)))))

pp_all_pri <- ppp(x = c(coords_pri[, 1], coords_ff[, 1]),
                  y = c(coords_pri[, 2], coords_ff[, 2]),
                  window = win,
                  marks = factor(c(rep("school", nrow(coords_pri)),
                                   rep("fastfood", nrow(coords_ff)))))

# 使用更細的距離序列（每10公尺）
r_vals <- seq(0, 3000, by = 10)

# 公立學校 → 速食店
F_pub <- envelope(pp_all_pub,
                  fun = Fest,
                  i = "school",
                  j = "fastfood",
                  nsim = 99,
                  r = r_vals,
                  verbose = FALSE)

# 私立學校 → 速食店
F_pri <- envelope(pp_all_pri,
                  fun = Fest,
                  i = "school",
                  j = "fastfood",
                  nsim = 99,
                  r = r_vals,
                  verbose = FALSE)


# 公立學校
plot(F_pub, main = "公立學校 → 速食店 的 F-cross 函數", legendargs = list(cex = 0.8))
abline(h = F_pub$theo, col = "blue", lty = 2)

# 私立學校
plot(F_pri, main = "私立學校 → 速食店 的 F-cross 函數", legendargs = list(cex = 0.8))
abline(h = F_pri$theo, col = "blue", lty = 2)

r_target <- 1000  # 想比較的距離

idx_pub <- which.min(abs(F_pub$r - r_target))
idx_pri <- which.min(abs(F_pri$r - r_target))

# 取得各自的 F_cross(1000)
F_pub_obs <- F_pub$obs[idx_pub]
F_pri_obs <- F_pri$obs[idx_pri]

cat("F(1000) 公立學校 → 速食店 =", F_pub_obs, "\n")

F(1000) 公立學校 → 速食店 = 0.4083 

cat("F(1000) 私立學校 → 速食店 =", F_pri_obs, "\n")

F(1000) 私立學校 → 速食店 = 0.5777631 

plot(F_pub$r, F_pub$obs, type = "l", col = "blue", lwd = 2,
     ylim = c(0, max(F_pub$obs, F_pri$obs)),
     xlab = "距離 r", ylab = "F_cross(r)",
     main = "F-cross 比較：公立 vs 私立")
lines(F_pri$r, F_pri$obs, col = "red", lwd = 2)

legend("bottomright", legend = c("公立學校", "私立學校"), col = c("blue", "red"), 
       lty = 1, lwd = 2)
abline(a = 0, b = 1 / max(F_pub$r), col = "gray", lty = 2)  # approx theo line

Part 3

方法一：倡導學校引導學生選擇健康餐盒

healthy_sf <- health_sf %>%
  mutate(行政區 = str_extract(Description, "行政區: [^<]+") %>% 
                   str_remove("行政區: "))

health_sf <- healthy_sf %>%
  filter(行政區 == "信義區")

schools_sf = st_join(SCHOOL_sf, Taipei_Vill_sf, join = st_within)%>%
  filter(TOWN %in% c("信義區"))

win_polygon <- st_as_sfc(st_bbox(health_sf))
win <- as.owin(win_polygon)

coords <- st_coordinates(health_sf)
healthyfood_ppp <- ppp(x = coords[, 1], y = coords[, 2], window = win)

r_seq <- seq(0, 5000, by = 100)

k_result <- Kest(healthyfood_ppp, r = r_seq)
k_3000 <- k_result$iso[which(k_result$r == 3000)]
l_3000 <- sqrt(k_3000 / pi) -3000

cat("K(3000) =", k_3000, "\n")

K(3000) = NA 

cat("L(3000) =", l_3000, "\n")

L(3000) = NA 

env <- envelope(healthyfood_ppp,
                fun = Kest,
                nsim = 99,
                r = r_seq,
                correction = "iso",
                savefuns = TRUE,
                savepatterns = TRUE,
                verbose = FALSE)

r_index <- which(env$r == 3000)
upper_L <- sqrt(env$hi[r_index] / pi) -3000
lower_L <- sqrt(env$lo[r_index] / pi) -3000

cat("95% CI for L(3000): [", lower_L, ",", upper_L, "]\n")

95% CI for L(3000): [ NA , NA ]

plot(env, . ~ r, main = "Ripley's k-function with 95% CI")
abline(v = 3000, col = "red", lty = 2)

coords_pub <- st_coordinates(schools_sf)
coords_ff  <- st_coordinates(health_sf)

win <- as.owin(st_as_sfc(st_bbox(schools_sf)))

pp_all_pub <- ppp(x = c(coords_pub[, 1], coords_ff[, 1]),
                  y = c(coords_pub[, 2], coords_ff[, 2]),
                  window = win,
                  marks = factor(c(rep("school", nrow(coords_pub)),
                                   rep("healthyfood", nrow(coords_ff)))))

r_vals <- seq(0, 3000, by = 2)

F_pub <- envelope(pp_all_pub,
                  fun = Fest,
                  i = "school",
                  j = "healthyfood",
                  nsim = 99,
                  r = r_vals,
                  verbose = FALSE)

plot(F_pub, main = "學校 → 健康餐盒店 的 F-cross 函數", legendargs = list(cex = 0.8))
abline(h = F_pub$theo, col = "blue", lty = 2)

方法二：補助學校周邊健康餐盒業者

#讀取資料
villages <- st_read("./data/Taipei_Vill.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Taipei_Vill' from data source `D:\R\spatial analysis\Final\data\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

schools <- st_read("./data/SCHOOL.shp", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `SCHOOL' from data source `D:\R\spatial analysis\Final\data\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

healthy_food <- st_read("./data/臺北市健康餐飲.kml", options="ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `臺北市健康餐飲' from data source 
  `D:\R\spatial analysis\Final\data\臺北市健康餐飲.kml' using driver `KML'
Simple feature collection with 234 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: 121.4936 ymin: 24.96722 xmax: 121.6181 ymax: 25.13751
Geodetic CRS:  WGS 84

#確認CRS
healthy_food <- st_transform(healthy_food, crs = st_crs(schools))

#建立 1500 公尺的緩衝區（以學校為中心）
schools_buffer <- st_buffer(schools, dist = 1500)

#找出每所學校 buffer 內的健康餐飲
school_food_join <- st_join(healthy_food, schools_buffer, join = st_within)


#統計每所學校有多少個健康餐飲據點
result <- school_food_join %>%
  st_drop_geometry() %>%
  group_by(Name.y) %>%
  summarise(HealthyFood_Count = n())

names(schools)

[1] "Type"     "Name"     "SID"      "index"    "geometry"

# [1] "編號" "NAME" "地址" "備註" "geometry"

names(result)

[1] "Name.y"            "HealthyFood_Count"

# [1] "Name.y" "HealthyFood_Count"

schools_summary <- left_join(schools, result, by = c("Name" = "Name.y"))

library(tmap)
base_lyr <- tm_shape(Taipei_Vill_sf) +
  tm_polygons(fill = "white", alpha = 1, border.col = "gray60") +
  tm_layout(frame = FALSE)

base_lyr+
tm_shape(schools_summary) +
  tm_dots(size = 0.5, col = "HealthyFood_Count", palette = "Greens", title = "1500m 內健康餐飲數") +
  tm_layout(title = "台北市學校周邊健康飲食分布")

villages <- st_transform(villages, crs = st_crs(schools))
schools_with_town <- st_join(schools, villages, join = st_within)
school_xinyi <- schools_with_town %>%
  filter(TOWN %in% c("信義區"))
school_xinyi_buffer <- st_buffer(school_xinyi, dist = 1500)

#對應健康餐飲點位
xinyi_school_food_join <- st_join(healthy_food, school_xinyi_buffer, join = st_within)

#統計每所學校的健康餐飲數量
xinyi_result <- xinyi_school_food_join %>%
  st_drop_geometry() %>%
  group_by(Name.y) %>%
  summarise(HealthyFood_Count = n())

#加回geometry
xinyi_summary <- left_join(school_xinyi, xinyi_result, by = c("Name" = "Name.y"))
base_lyr <- tm_shape(xinyi_sf) +
  tm_polygons(fill = "white", alpha = 1, border.col = "gray60") +
  tm_layout(frame = FALSE)

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

base_lyr+
tm_shape(xinyi_summary) +
  tm_dots(size = 0.5, col = "HealthyFood_Count", palette = "Greens", title = "1500m 內健康餐飲數") +
  tm_layout(title = "信義區學校周邊健康飲食分布")

[tm_dots()] Argument `title` unknown.[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run `cols4all::c4a_gui()` to explore
them. The old palette name "Greens" is named "brewer.greens"Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

library(tmap)

library(sf)
library(dplyr)
library(tmap)



# 投影轉換：確保所有資料都是一樣的 CRS
healthy_food <- st_transform(healthy_food, st_crs(schools))
villages <- st_transform(villages, st_crs(schools))

# 🔧 處理 MULTIPOINT + Z 維度問題
healthy_food <- healthy_food %>%
  st_cast("POINT") %>%
  st_zm(drop = TRUE)
# 找出信義區學校
schools_with_town <- st_join(schools, villages, join = st_within)
school_xinyi <- schools_with_town %>%
  filter(TOWN %in% c("信義區"))

# 建立緩衝區
school_xinyi_buffer <- st_buffer(school_xinyi, dist = 500)

# 加入健康餐飲點位（找落在 buffer 裡的）
xinyi_school_food_join <- st_join(healthy_food, school_xinyi_buffer, join = st_within)

# 統計每間學校的餐飲點數量
xinyi_result <- xinyi_school_food_join %>%
  st_drop_geometry() %>%
  group_by(Name.y) %>%
  summarise(HealthyFood_Count = n())

#加geometry
xinyi_summary <- left_join(school_xinyi, xinyi_result, by = c("Name" = "Name.y"))

#視覺化
tmap_mode("plot") 

ℹ tmap mode set to "plot".

# 地圖主體與圖層
base_lyr <- tm_shape(xinyi_sf) +
  tm_polygons(fill = "white", alpha = 1, border.col = "gray60") +
  tm_layout(frame = FALSE)

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

base_lyr+
tm_shape(school_xinyi_buffer) +
  tm_polygons(
    alpha = 0.2,
    border.col = "blue",
    border.lwd = 1
  ) +
tm_shape(healthy_food) +
  tm_dots(
    size = 0.2,
    fill = "red",
    fill.alpha = 0.7,
    legend.show = TRUE
  ) +
tm_shape(school_xinyi) +
  tm_dots(
    size = 0.06,
    fill = "black",
    col = "white"
  ) +
  tm_text("Name", size = 0.6, col = "black", just = "top") +
tm_title("信義區學校及其 1500 公尺內健康餐飲分布") +
tm_layout(legend.outside = TRUE)

[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.[tm_polygons()] Argument `border.lwd` unknown.[tm_dots()] Arguments `fill.alpha` and `legend.show` unknown.[v3->v4] `tm_text()`: migrate the layer options 'just' to 'options = opt_tm_text(<HERE>)'

#install.packages("gt")
library(gt)

警告： 套件 ‘gt’ 是用 R 版本 4.4.3 來建造的
載入套件：‘gt’

下列物件被遮斷自 ‘package:Hmisc’:

    html

下列物件被遮斷自 ‘package:tmap’:

    metro

library(tmap)
library(gt)

xinyi_result %>%
  gt() %>%
  tab_header(
    title = "信義區學校周邊健康飲食點統計",
    subtitle = "500 公尺範圍內的健康餐飲據點數"
  ) %>%
  cols_label(
    Name.y = "學校名稱",
    HealthyFood_Count = "健康餐飲數量"
  ) %>%
  fmt_number(
    columns = HealthyFood_Count,
    decimals = 0
  )

信義區學校周邊健康飲食點統計
500 公尺範圍內的健康餐飲據點數
學校名稱	健康餐飲數量
三興國小	3
信義國小	1
光復國小	6
博愛國小	1
吳興國小	1
永吉國小	3
福德國小	1
興雅國小	1
NA	217

tmap_mode("view") 

ℹ tmap mode set to "view".

tm_shape(school_xinyi_buffer) +
  tm_polygons(
    alpha = 0.2,
    border.col = "blue",
    border.lwd = 1
  ) +
tm_shape(healthy_food) +
  tm_dots(
    size = 0.2,
    fill = "red",
    fill.alpha = 0.7,
    legend.show = TRUE
  ) +
tm_shape(school_xinyi) +
  tm_dots(
    size = 0.06,
    fill = "black",
    col = "white"
  ) +
  tm_text("Name", size = 0.6, col = "black", just = "top") +
tm_title("信義區學校及其 1500 公尺內健康餐飲分布") +
tm_layout(legend.outside = TRUE)

── tmap v3 code detected ──────────────────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.[tm_polygons()] Argument `border.lwd` unknown.[tm_dots()] Arguments `fill.alpha` and `legend.show` unknown.[v3->v4] `tm_text()`: migrate the layer options 'just' to 'options = opt_tm_text(<HERE>)'

方法三：綠地派對

• 找離central feature最近的公園

find_CentralFeature = function(xy_data_frame){
  CF = calc_cf(id=1, points = xy_data_frame)
  CF.x = CF$ATTRIBUTES$CF.x
  CF.y = CF$ATTRIBUTES$CF.y
  CF.coor = c(CF.x, CF.y)
  CF_sf = CF.coor %>% st_point %>% st_sfc %>% st_sf
  return(CF_sf)
}

xinyi_index =  Taipei_Vill_sf$TOWN == "信義區"
xinyi_sf = Taipei_Vill_sf[xinyi_index,] 

Park_sf = st_read("./data/L0101-5.geojson") 

Reading layer `L0101-5' from data source `D:\R\spatial analysis\Final\data\L0101-5.geojson' using driver `GeoJSON'
Simple feature collection with 880 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 296575.9 ymin: 2763553 xmax: 313042.4 ymax: 2782160
Projected CRS: TWD97 / TM2 zone 121

Park_xin_poly_sf = st_join(Park_sf, Taipei_Vill_sf, join = st_within)%>%
  filter(TOWN %in% c("信義區"))
Park_xinyi_sf = st_centroid(Park_sf) %>% 
  st_join(Taipei_Vill_sf, join = st_within)%>%
  filter(TOWN %in% c("信義區"))

school_xinyi_sf = st_join(SCHOOL_sf, Taipei_Vill_sf, join = st_within)%>%
  filter(TOWN %in% c("信義區"))

school_coor = st_coordinates(school_xinyi_sf)
school_xinyi_sf$x = school_coor[, 1]
school_xinyi_sf$y = school_coor[, 2]
school_xinyi_coor_df = as.data.frame(school_xinyi_sf[,14:15])
CF_sf = find_CentralFeature(school_xinyi_coor_df[,1:2])
st_crs(CF_sf) = st_crs(SCHOOL_sf)

distances = st_distance(CF_sf, Park_xinyi_sf)
closest_index = which.min(distances)
closest_park = Park_xinyi_sf[closest_index, ]
closest_park_poly <- st_join(closest_park, Park_xin_poly_sf, join = st_intersects)

• 畫地圖

library(tmap)
tmap_mode("plot")

ℹ tmap mode set to "plot".

# 地圖主體與圖層
base_lyr <- tm_shape(xinyi_sf) +
  tm_polygons(fill = "white", alpha = 1, border.col = "gray60") +
  tm_layout(frame = FALSE)

── tmap v3 code detected ──────────────────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

# 學校點位
school_lyr <- tm_shape(school_xinyi_sf) +
  tm_dots(fill = "blue", size = 0.2, fill.alpha = 0.8,
          title = "學校")

[tm_dots()] Arguments `fill.alpha` and `title` unknown.

# 中央點 CF
cf_lyr <- tm_shape(CF_sf) +
  tm_symbols(shape = 21, fill = "red", size = 0.6,
             title = "學校重心點")

[tm_symbols()] Argument `title` unknown.

# 所有公園（背景）
park_lyr <- tm_shape(Park_xin_poly_sf) +
  tm_polygons(fill = "lightgreen", alpha = 0.4, border.col = NA,
              legend.show = FALSE)  # 不讓背景公園出現在圖例

[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.[v3->v4] `tm_polygons()`: use `fill.legend = tm_legend_hide()` instead of `legend.show = FALSE`.

# 最近的那一塊公園（重點標記）
closest_park_poly <- st_join(closest_park, Park_xin_poly_sf, join = st_nearest_feature)

closest_park_lyr <- tm_shape(closest_park_poly) +
  tm_polygons(fill = "forestgreen", alpha = 0.7,
              border.col = "darkgreen", lwd = 2,
              title = "最近的公園")

[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the visual variable `fill` namely
'title' to 'fill.legend = tm_legend(<HERE>)'

# 組合地圖並加上標題與圖例設定

final_map <- base_lyr + park_lyr + closest_park_lyr + school_lyr + cf_lyr +
  tm_add_legend(type = "fill", labels = c("信義區公園", "綠地派對公園"), 
                col = c("lightgreen", "forestgreen")) +
  tm_add_legend(type = "symbol", labels = c("學校Central Feature"), col = "red", shape = 21, size = 0.6) +
  tm_add_legend(type = "symbol", labels = c("學校"), col = "blue", shape = 21, size = 0.2) +
  tm_scalebar(position = c("left", "bottom")) + 
  tm_compass(position = c("right", "top")) +
  tm_layout(
    main.title = "信義區：綠地派對舉辦地點",
    main.title.size = 1.2,
    legend.position = c("left", "bottom"),
    legend.bg.color = "white",
    legend.bg.alpha = 0.7,
    frame = FALSE
  )

[v3->v4] `tm_add_legend()`: use `type = "polygons"` instead of `type = "fill"`.[v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of polygons.[v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type = "symbol"`.[v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of symbols.[v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type = "symbol"`.[v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of symbols.[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

print(final_map)
tmap_save(final_map, filename = "xinyi_map.jpg", width = 2000, height = 1600, dpi = 300)