G12

rm(list = ls())

setwd("C:/Users/618am/OneDrive/桌面/空間分析期末報告")

library(sf)
library(tmap)
library(ggplot2)
library(dplyr)
library(units)
library(reshape2)
library(aspace)
library(spatstat)
library(spdep)
library(tidyverse)
library(pals)
library(cartography)
library(readr)

Popn_TWN2 <- st_read("C:/Users/618am/OneDrive/桌面/空間分析期末報告/Popn_TWN2", options = "ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Popn_TWN2' from data source 
  `C:\Users\618am\OneDrive\桌面\空間分析期末報告\Popn_TWN2' using driver `ESRI Shapefile'
Simple feature collection with 368 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -197572.5 ymin: 2295201 xmax: 606875.6 ymax: 2919551
Projected CRS: TWD97 / TM2 zone 121

SCHOOL <- st_read("C:/Users/618am/OneDrive/桌面/空間分析期末報告/SCHOOL", options = "ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `SCHOOL' from data source 
  `C:\Users\618am\OneDrive\桌面\空間分析期末報告\SCHOOL' 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

Tpe_Fastfood <- st_read("C:/users/618am/OneDrive/桌面/空間分析期末報告/Tpe_Fastfood", options = "ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Tpe_Fastfood' from data source 
  `C:\Users\618am\OneDrive\桌面\空間分析期末報告\Tpe_Fastfood' 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

TPE_Vill <- st_read("C:/Users/618am/OneDrive/桌面/空間分析期末報告/Taipei_Vill/Taipei_Vill.shp", options = "ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Taipei_Vill' from data source 
  `C:\Users\618am\OneDrive\桌面\空間分析期末報告\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

TPE_Village_basemap <- tm_shape(TPE_Vill) +
  tm_polygons(col = "grey") +
  tm_layout(title = "臺北市公私立國小與速食店點位分布圖")

[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

SCHOOL_map <- tm_shape(SCHOOL) +
  tm_dots(fill = "blue", size = 0.25)

Fastfood_map <- tm_shape(Tpe_Fastfood) +
  tm_dots(fill = "red", size = 0.25)

TPE_Village_basemap + SCHOOL_map + Fastfood_map + tm_add_legend(type = "symbol",
              labels = c("學校", "速食店"),
              col = c("blue", "red"),
              title = "類別") + tm_scale_bar() + tm_compass() + tm_layout()

[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.! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.

Part 1: 不同行政區的速食店有不同的地理分布特徵嗎？

問答題

圖 1 是某城市各村里速食店密度的 Moran Scatter Plot。

(1)

簡述該圖表的 X 軸與 Y 軸所表示的意涵，以及其空間型態。

X 軸代表各空間單位的速食店密度，經過標準化處理。數值越高，表示該區域速食店密度高於平均。 Y 軸表示「鄰近區域的速食店密度」的平均值，也經過標準化。這是計算空間自相關的關鍵：看某一區的速食店密度與其鄰近區域的密度是否一致。

(2)

請簡述該圖趨勢線(Fit)的斜率值 (slope =0.982)，所代表的意涵。

數值代表空間自相關的強度， 趨勢線的斜率近似 1，表示變數（速食店密度）與其鄰近區域的值之間有非常強的正向線性關係。

(3)

關於該城市各村里速食店密度的 Local Moran’s I 統計量，請問其統計量可能會大於 0、等於 0 或小於 0，並說明判斷的理由。

Local Moran’s I 可能大於 0（高-高、低-低聚集）、小於 0（離群區，如高-低、低-高），或等於 0（無空間關聯）。即使整體呈現正向空間自相關，局部仍可能存在異常點與離群現象，因此需使用 Local Moran’s I 進行細部判斷。

實作題

(1)

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

# 定義 A 區、B 區行政區名稱
A_towns <- c("文山區","大安區","中正區")
B_towns <- c("信義區","南港區","松山區")

# 1. 建立 1km 緩衝區 around each fastfood 
fast_buf <- st_buffer(Tpe_Fastfood, dist = set_units(1, km))

# 2. 空間連結：每家速食店緩衝內的學校數量
ff_school <- st_join(fast_buf, SCHOOL, join = st_contains)
counts <- ff_school %>% 
  st_drop_geometry() %>% 
  group_by(ID) %>%               
  summarise(school_n = n())

# 3. 對速食店附加其所屬行政區
ff_town <- st_join(Tpe_Fastfood, TPE_Vill, join = st_within)
ff_stats <- left_join(ff_town, counts, by = "ID")

# 4. 加入 A/B 類別標籤
ff_stats <- ff_stats %>%
  mutate(area_group = case_when(
    TOWN %in% A_towns ~ "A區",
    TOWN %in% B_towns ~ "B區",
    TRUE              ~ NA_character_
  )) %>%
  filter(!is.na(area_group))

# 5. 列出 A、B 區中每家速食店 的 school_n 向量，
#     並做 t.test（或 Wilcox test）
A_vals <- ff_stats$school_n[ff_stats$area_group=="A區"]
B_vals <- ff_stats$school_n[ff_stats$area_group=="B區"]

# 虛無假設 H0: μ_A = μ_B；對立假設 H1: μ_A ≠ μ_B
t_res <- t.test(A_vals, B_vals, var.equal = FALSE)
t_res

    Welch Two Sample t-test

data:  A_vals and B_vals
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 

虛無假設（H₀）：A 區與 B 區的速食店服務範圍內涵蓋的學校數平均值相同（μ₁ = μ₂）。 對立假設（H₁）：A 區與 B 區的速食店涵蓋學校數平均值不同（μ₁ ≠ μ₂）。

A區與B區速食店涵蓋的學校數平均值無顯著差異（p = 0.1163 > 0.05），無法拒絕虛無假設。顯示兩區速食店服務效能在統計上無明顯差異。

(2)

以 500 公尺方格的空間單位建立網格，依照 contiguity 鄰近定義，計算並繪製速食店總數的 Moran’s correlogram，解讀其圖表資訊，並評估速食店的空間影響範圍。

# 1. 建立 500 公尺格網，並裁切至台北市行政區內部
grid_500 <- st_make_grid(TPE_Vill, cellsize = 500, what = "polygons") %>%
  st_sf() %>%
  st_intersection(st_union(TPE_Vill))

grid_500$grid_id <- 1:nrow(grid_500)

# 2. 空間連結速食店到格網，統計每格速食店數量
grid_join <- st_join(Tpe_Fastfood, grid_500, join = st_within)
count_grid <- grid_join %>%
  group_by(grid_id) %>%
  summarise(n = n())

# 3. 將計數結果合併回格網資料，將 NA 補為 0（無速食店）
count_grid_nogeo <- count_grid %>% st_drop_geometry()
grid_500 <- left_join(grid_500, count_grid_nogeo, by = "grid_id") %>%
  mutate(n = replace_na(n, 0))

# 4. 建立格網之鄰接矩陣（以 Queen 鄰接為準）
nb_grid <- poly2nb(grid_500, queen = TRUE)
lw_grid <- nb2listw(nb_grid, style = "W", zero.policy = TRUE)

# 5. 提取速食店數值向量
n_vector <- grid_500$n

# 6. 計算 Moran’s I Correlogram（鄰近階數從 1~6）
moran_correlogram <- sp.correlogram(
  nb_grid, 
  n_vector, 
  order = 6, 
  method = "I", 
  style = "W", 
  zero.policy = TRUE
)

# 7. 繪製 correlogram 結果圖
plot(moran_correlogram, main = "Moran's I Correlogram of Fast Food Counts (500m Grid)")

Moran’s I 在第一階層鄰近（Lag 1）顯著為正，顯示速食店有明確的局部群聚性，但隨距離增加自相關性下降，影響範圍約在 500–1000 公尺內。

tm_shape(grid_500) +
  tm_polygons(col = "n", palette = "YlOrRd", title = "速食店數") +
  tm_layout(title = "臺北市 500m 網格中速食店分布")

── 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(title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run `cols4all::c4a_gui()` to
explore them. The old palette name "YlOrRd" is named "brewer.yl_or_rd"Multiple palettes called "yl_or_rd" found: "brewer.yl_or_rd", "matplotlib.yl_or_rd". The first one, "brewer.yl_or_rd", is returned.

(3)

依照前一小題的網格為空間單位，依照 contiguity 鄰近定義，計算 Gi*統計量，分別繪製速食店家數與學校在 0.05 顯著水準的熱區分佈（需進行多重檢定的修正）。

school_join <- st_join(SCHOOL, grid_500, join = st_within)
school_count <- school_join %>%
  group_by(grid_id) %>%
  summarise(schools = n())

school_count_nogeo <- school_count %>% st_drop_geometry()
grid_500 <- left_join(grid_500, school_count_nogeo, by = "grid_id") %>%
  mutate(schools = replace_na(schools, 0))

nb_grid <- poly2nb(grid_500, queen = TRUE)
lw_grid <- nb2listw(nb_grid, style = "W", zero.policy = TRUE)

# 計算 Gi* 統計量
Gi_fast <- localG(grid_500$n, lw_grid)
Gi_school <- localG(grid_500$schools, lw_grid)

# 計算 p-value 並進行 FDR 修正（雙尾）
grid_500$Gi_fast <- as.numeric(Gi_fast)
grid_500$p_fast <- 2 * pnorm(-abs(grid_500$Gi_fast))
grid_500$FDR_fast <- p.adjust(grid_500$p_fast, method = "fdr")
grid_500$Gi_school <- as.numeric(Gi_school)
grid_500$p_school <- 2 * pnorm(-abs(grid_500$Gi_school))
grid_500$FDR_school <- p.adjust(grid_500$p_school, method = "fdr")

# 熱區分類標籤（顯著性標示）
grid_500 <- grid_500 %>%
  mutate(hotspot_fast = case_when(
    FDR_fast < 0.05 & Gi_fast > 0 ~ "High-High",
    FDR_fast < 0.05 & Gi_fast < 0 ~ "Low-Low",
    TRUE ~ "Not significant"
  )) %>%
  mutate(hotspot_school = case_when(
    FDR_school < 0.05 & Gi_school > 0 ~ "High-High",
    FDR_school < 0.05 & Gi_school < 0 ~ "Low-Low",
    TRUE ~ "Not significant"
  ))

# 熱區圖繪製（速食店）
tm_shape(grid_500) +
  tm_fill("hotspot_fast",
          palette = c("red", "white", "gray80"),
          title = "速食店熱區",
          labels = c("High-High", "Low-Low", "Not significant")) +
  tm_borders() +
  tm_layout(title = "速食店 Gi* 熱區分析 (FDR < 0.05)")

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the visual variable
`fill` namely 'palette' (rename to 'values'), 'labels' to fill.scale = tm_scale(<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(title = )`警告： labels do not have the same length as levels, so they are repeated

# 熱區圖繪製（學校）
tm_shape(grid_500) +
  tm_fill("hotspot_school",
          palette = c("red", "white", "gray80"),
          title = "學校熱區",
          labels = c("High-High", "Low-Low", "Not significant")) +
  tm_borders() +
  tm_layout(title = "學校 Gi* 熱區分析 (FDR < 0.05)")

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the visual variable
`fill` namely 'palette' (rename to 'values'), 'labels' to fill.scale = tm_scale(<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(title = )`警告： labels do not have the same length as levels, so they are repeated

速食店與學校的熱區皆集中於臺北市西南部，顯示兩者空間分布具有重疊群聚現象。該區可視為學生容易接觸速食的高風險環境。

Part 2: 學校附近的速食店數量真的比較多嗎?

(1)

請參考 Transportation Research Part A, 45 (2011) 640–652（HW-09 的研讀教材）所使用 F(d)函數的作法，比較分別位於大安區、文山區與信義區的學校，到鄰近速食店的空間特徵。

SCHOOL <- st_join(SCHOOL, TPE_Vill[, c("TOWN")], join = st_within)
table(SCHOOL$TOWN)

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

districts <- c("大安區", "文山區", "信義區")
schools_by_district <- SCHOOL %>% filter(TOWN %in% districts)

# 計算距離最近速食店的距離
dist_list <- lapply(districts, function(dist){
  sch <- schools_by_district %>% filter(TOWN == dist)
  m <- st_distance(sch, Tpe_Fastfood)
  dmin <- apply(m, 1, min) %>% set_units("m") %>% drop_units()
  data.frame(
    district = dist,
    dist_m   = dmin
  )
})
df_dist <- bind_rows(dist_list)

# 繪製 F(d) ECDF 圖
ggplot(df_dist, aes(x = dist_m, colour = district)) +
  stat_ecdf(size = 1.2) +
  scale_x_continuous(limits = c(0, 2000),
                     breaks = seq(0, 2000, 500),
                     labels = function(x) paste0(x, "m")) +
  labs(
    title = "臺北市三區學校到最近速食店的 F(d) 距離分布",
    x = "距離 d（公尺）",
    y = "F(d) = 學校距離 ≤ d 的比例",
    colour = "行政區"
  ) +
  theme_minimal() +
  theme(text = element_text(family = "Heiti TC"), legend.position = "bottom")

# KS檢定（兩兩組合）
ks_results <- combn(districts, 2, simplify = FALSE, FUN = function(pair){
  d1 <- df_dist %>% filter(district == pair[1]) %>% pull(dist_m)
  d2 <- df_dist %>% filter(district == pair[2]) %>% pull(dist_m)
  kt <- ks.test(d1, d2)
  data.frame(
    group1 = pair[1],
    group2 = pair[2],
    statistic = kt$statistic,
    p_value  = kt$p.value
  )
})
ks_df <- bind_rows(ks_results)
print(ks_df)

信義區學校距速食店最近，文山區最遠，顯示信義區速食店密度高、可及性強；文山區則相對分散，學生暴露速食環境風險較低。

(2)

利用 Ripley’s K function, k(d)計算速食店在 distance (d) = 3000 公尺時的 K(3000) 以及L(3000)，並使用 Monte Carlo significance test 計算 L(3000)的 95%信賴區間。

# 1. 準備 ppp 資料（速食店點與分析範圍）
# 轉換為 spatstat 分析視窗（以 TPE_Vill 範圍為主）
win <- as.owin(st_union(TPE_Vill))
fast_xy <- st_coordinates(Tpe_Fastfood)
ppp_fast <- ppp(x = fast_xy[,1], y = fast_xy[,2], window = win)

警告： data contain duplicated points

# 2. 計算 K(d) 與 L(d)（使用 border 修正）
K_result <- Kest(ppp_fast, correction = "border")
L_result <- Lest(ppp_fast, correction = "border")

# 3. 提取 d = 3000 時的 K 與 L 值
target_d <- 3000
k_index  <- which.min(abs(K_result$r - target_d)) 

k_r_3000 <- K_result$r[k_index]                    
k_3000   <- K_result$border[k_index]               
l_3000   <- L_result$border[k_index]               

cat("K(", k_r_3000, ") =", round(k_3000, 2), "\n")

K( 2998.733 ) = 67640147 

cat("L(", k_r_3000, ") =", round(l_3000, 2), "\n")

L( 2998.733 ) = 4640.1 

# 4. Monte Carlo 模擬信賴區間（99 次 CSR 模擬）
set.seed(1234)
L_env <- envelope(
  ppp_fast,
  fun = Lest,
  nsim = 99,
  correction = "border",
  global = FALSE,
  simulate = expression(rpoispp(lambda = intensity(ppp_fast), win = win))
)

Generating 99 simulations by evaluating expression  ...
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.

L_3000_obs <- L_env$obs[k_index]
L_3000_lo  <- L_env$lo[k_index]
L_3000_hi  <- L_env$hi[k_index]

cat("L(", k_r_3000, ") =", round(L_3000_obs, 2), 
    "\n95% CI = [", round(L_3000_lo, 2), ",", round(L_3000_hi, 2), "]\n")

L( 2998.733 ) = 4640.1 
95% CI = [ 2357.99 , 3418.9 ]

# 5. 繪圖：Ripley's L with 信賴區間與標示線 
plot(L_env,
     main = "Ripley's L function with 95% Monte Carlo Envelope",
     xlab = "距離 d (公尺)", ylab = "L(d)",
     xlim = c(0, 4000),
     legend = FALSE)
abline(h = 0, col = "gray30", lty = 2)             

abline(v = k_r_3000, col = "red", lty = 3)

# 在圖上加上實際值與信賴區間的點
points(x = k_r_3000, y = L_3000_obs, pch = 16, col = "black")

segments(x0 = k_r_3000, x1 = k_r_3000, y0 = L_3000_lo, y1 = L_3000_hi, col = "darkblue", lwd = 2)

L(3000) 值顯著高於 CSR 模擬信賴區間，表示在 3 公里尺度下速食店存在顯著空間群聚，顯示其分布集中且非隨機。

(3)

利用 Bivariate F function 分析方法，比較公立 vs.私立學校，到鄰近速食店的空間特徵，並說明哪一種類型學校的附近，會有較多的速食店？（顯著水準=0.05）

table(SCHOOL$Type)

公立 私立 
  51   97 

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

# 計算兩群學校到最近速食店的距離
get_nearest_dist <- function(school_sf, fast_sf){
  dist_matrix <- st_distance(school_sf, fast_sf)
  dmin <- apply(dist_matrix, 1, min) %>% set_units("m") %>% drop_units()
  return(dmin)
}

dist_pub <- get_nearest_dist(school_public, Tpe_Fastfood)
dist_pri <- get_nearest_dist(school_private, Tpe_Fastfood)

# 合併資料
df_dist <- data.frame(
  dist_m = c(dist_pub, dist_pri),
  Type   = c(rep("公立", length(dist_pub)), rep("私立", length(dist_pri)))
)

# 繪圖：F(d) 比較（使用 ECDF）
ggplot(df_dist, aes(x = dist_m, colour = Type)) +
  stat_ecdf(size = 1.2) +
  scale_x_continuous(limits = c(0, 2000),
                     breaks = seq(0,2000,500),
                     labels = function(x) paste0(x, "m")) +
  labs(
    title = "公私立學校到速食店的最近距離 F(d) 比較",
    x     = "距離 d（公尺）",
    y     = "F(d) = 學校距離 ≤ d 的比例",
    colour= "學校類型"
  ) +
  theme_minimal() +
  theme(text = element_text(family = "Heiti TC"), legend.position = "bottom")

# KS 檢定：比較距離分布是否顯著不同
ks_test <- ks.test(dist_pub, dist_pri)
print(ks_test)

    Exact two-sample Kolmogorov-Smirnov test

data:  dist_pub and dist_pri
D = 0.35052, p-value = 0.0003497
alternative hypothesis: two-sided

公立學校 F(d)曲線上升較快，表示距離速食店較近、可及性更高。私立學校相對距離較遠，顯示公立學校周邊速食密度較高。

Part 3: 校園健康生活計畫

根據臺北市衛生局於113年公布之健康元素餐飲業者清冊，記錄每間健康飲食餐廳點位，並針對臺北市每間國小方圓一公里進行環域分析，記錄每間學校環域內健康飲食店家，提供學童在外食上更多元且營養豐富的消費選擇，並針對環域內未記錄到健康飲食店家之學校制定進一步校園健康飲食計畫。

healthy_restaurant <- read.csv("C:/Users/618am/OneDrive/桌面/空間分析期末報告/臺北市健康飲食餐廳資料表 - healthy_food_shops_with_address.csv", fileEncoding = "UTF-8")
head(healthy_restaurant)

healthy_geom <- st_as_sf(
  healthy_restaurant,
  coords = c("經度", "緯度"),
  crs = 4326
)

healthy_geom_3826 <- st_transform(healthy_geom, 3826)
SCHOOL_buffer <- st_buffer(SCHOOL, dist = 1000)

school_healthy <- st_join(SCHOOL_buffer, healthy_geom_3826, join = st_intersects, left = TRUE)
summary_table <- school_healthy %>%
  st_drop_geometry() %>%
  group_by(Name) %>%   
  summarise(店家數 = sum(!is.na(商家名稱)))  

school_healthy_info <- school_healthy %>%
  st_drop_geometry()
SCHOOL_enriched <- SCHOOL %>%
  left_join(school_healthy_info, by = "SID")

healthy_geom_df <- healthy_geom_3826 %>%
  mutate(healthy_geom = geometry) %>%
  st_drop_geometry() %>%
  select(商家名稱, healthy_geom)

school_healthy_distance <- SCHOOL_enriched %>%
  left_join(healthy_geom_df, by = "商家名稱") %>%
  mutate(距離 = as.numeric(st_distance(geometry, healthy_geom, by_element = TRUE)))

schools_no_restaurant <- school_healthy %>% filter(is.na(商家名稱))
schools_no_restaurant$Name

 [1] "民權國小" "三民國小" "大直國小" "永安國小" "濱江國小" "博嘉國小" "萬芳國小" "力行國小"
 [9] "東湖國小" "明湖國小" "南湖國小" "麗湖國小" "陽明國小" "社子國小" "富安國小" "劍潭國小"
[17] "溪山國小" "平等國小" "百齡國小" "雙溪國小" "葫蘆國小" "湖田國小" "泉源國小" "大屯國小"
[25] "湖山國小"

一公里環域內無健康飲食餐廳學校

TPE_Village_basemap_2 <- tm_shape(TPE_Vill) +
  tm_polygons(col = "grey") +
  tm_layout(title = "臺北市公私立國小與健康飲食餐廳與速食店點位分布圖")

[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

SCHOOL_map <- tm_shape(SCHOOL) +
  tm_dots(fill = "blue", size = 0.25)

Healthy_food_map <- tm_shape(healthy_geom_3826) +
  tm_dots(fill = "green", size = 0.25)

TPE_Village_basemap_2 + SCHOOL_map + Healthy_food_map+ Fastfood_map + 
  tm_add_legend(type = "symbol",
              labels = c("學校", "健康飲食餐廳", "速食店"),
              col = c("blue", "green", "red"),
              title = "類別") + 
  tm_scale_bar() + 
  tm_compass() + 
  tm_layout()

[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.! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.

tm_shape(TPE_Vill) + tm_polygons(col = "grey") + 
  SCHOOL_map + Healthy_food_map + 
  tm_shape(SCHOOL_buffer) + 
  tm_borders(col = "black", lwd = 1) +
  tm_add_legend(type = "symbol", labels = c("學校", "健康飲食餐廳"), col = c("blue", "green"), title = "類別")

[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.

school_healthy_df <- st_drop_geometry(school_healthy_distance)
write_csv(school_healthy_df, "school_healthy_result_with_distance.csv")

學校與速食店熱區重疊篩選及健康餐車推薦設置地點

grid_500 <- grid_500 %>%
  mutate(
    joint = case_when(hotspot_fast == "High-High" & hotspot_school == "High-High" ~ "Both High",
                      hotspot_fast == "High-High" ~ "Fast-food Only",
                      hotspot_school == "High-High" ~ "Schools Only",
                      TRUE ~ "Neither"),
    joint = factor(joint, levels = c("Both High","Fast-food Only", "Schools Only", "Neither"))
  )
buf300 <- st_buffer(SCHOOL, dist = 300)
candidates <- grid_500 %>%
  filter(joint == "Both High") %>%
  filter(lengths(st_intersects(.,buf300))>0)
healthy_trucks_for_both_high <- st_centroid(st_geometry(candidates)) %>% st_as_sf()
healthy_trucks_for_no_restaurant <- st_centroid(st_geometry(schools_no_restaurant)) %>% st_as_sf()

tm_shape(grid_500) + tm_polygons(col = "joint", fill.scale = tm_scale(values = c("red","orange","blue","gray80")), fill.legend = tm_legend(title = "Hotspot Category"),
                                 border.col = "white",
                                 col_alpha = 0.3) +
  tm_layout(main.title = "School vs Fast-Food Bivariate Hotspots",
            legend.outside = TRUE,
            frame = FALSE)

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

尋找位於高風險區域學校

candidates_buf <- st_buffer(candidates, dist = 300)
school_near <- SCHOOL %>% filter(lengths(st_intersects(., candidates_buf))>0)
school_names <- unique(school_near$Name)
school_names

[1] "信義國小" "仁愛國小" "南門國小" "北教實小" "蓬萊國小" "日新國小" "雙蓮國小" "大同國小"

tm_shape(TPE_Vill) +
  tm_borders(col = "gray70") +

tm_shape(candidates) +
  tm_polygons(
    col         = "red",
    fill_alpha  = 0.2,
    border.col  = "red",
    fill.legend = "High-risk Grid (Both Hotspots)"
  ) +

tm_shape(healthy_trucks_for_no_restaurant) +
  tm_symbols(
    col         = "green",
    shape       = 21,
    size        = 0.25,
    border.col  = "darkgreen",
    legend.show = TRUE,
    fill.legend = "Recommended Food Truck Site"
  ) +
tm_add_legend(type = "fill", 
              labels = "High-risk Grid", 
              col = "red", 
              alpha = 0.2) +

tm_add_legend(type = "symbol", 
              labels = "Recommended Food Truck Site", 
              col = "green", 
              shape = 21) +
tm_layout(
  main.title      = "Recommended Healthy Food Truck Sites",
  main.title.size = 1.2,
  main.title.position = "center",
  legend.outside  = TRUE,
  legend.title.size = 1.0,
  legend.text.size  = 0.8,
  frame            = FALSE
)

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────
[tm_symbols()] Argument `legend.show` unknown.[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_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

針對高風險地區及方圓一公里範圍內無健康飲食餐廳之學校，提供健康飲食餐車建議設置地點。