台北市國小校園周邊健康生活計畫

rm(list = ls())
setwd("~/Documents/113-2/空間分析/final")

library(sf);library(tidyverse);library(readxl);library(jsonlite);library(spatstat);library(ggplot2);library(ggspatial);library(showtext);library(sfdep);library(tmap);library(dplyr)

# 設定台北市邊界資料（若有 shp）
taipei <- st_read('./data/Taipei_village/Taipei_village.shp', options = 'ENCODING = BIG5', quiet = T)
taipei <- st_transform(taipei, 4326)
taipei <- taipei %>% group_by(TOWNNAME) %>% summarise()

# 學校資料（自行提供）
School <- st_read("./data/SCHOOL/SCHOOL.shp", quiet = TRUE) # 含公私立屬性欄位
School <- st_transform(School, 4326)

# 速食店資料
fastfood <- st_read("./data/Tpe_Fastfood/Tpe_Fastfood.shp", quiet = TRUE)
fastfood <- st_transform(fastfood, 4326)

# 公園資料
Park <- fromJSON("./data/臺北市公園基本資料.json") %>%
  st_as_sf(coords = c("pm_Longitude", "pm_Latitude"), crs = 4326)

# 醫院資料
Hospital <- read_excel("./data/台北市醫院清冊1130523(含經緯度).xlsx") %>%
  st_as_sf(coords = c("經度", "緯度"), crs = 4326)

# 台北診所資料
clinic <- st_read("./data/台北診所資料.geojson", quiet = TRUE)

# 運動中心
Sports_Center <- read.csv("./data/臺北市運動中心.csv", encoding = "UTF-8") %>%
  st_as_sf(coords = c("經度", "緯度"), crs = 4326)

# 健康飲食地點
health_food <- st_read("./data/health/health.shp", quiet = TRUE)
health_food <- st_transform(health_food, 4326)

# 河川資料
river <- st_read("./data/RIVERPOLY/riverpoly/riverpoly.shp", quiet = TRUE)
river <- st_transform(river, 4326)
river <- river[river$RIVER_NAME == "淡水河" | river$RIVER_FROM == "淡水河",]
river <- river[!is.na(river$RIVER_NAME), ]

# 登山步道資料
trail_start <- read.csv("./data/登山步道.csv", encoding = "UTF-8") 
trail_start <- trail_start[-32, ]
trail_start <- st_as_sf(trail_start, coords = c("起點經度座標", "起點緯度座標"), crs = 4326)

trail_end <- read.csv("./data/登山步道.csv", encoding = "UTF-8") 
trail_end <- trail_end[-32, ]
trail_end <- st_as_sf(trail_start, coords = c("終點經度座標", "終點緯度座標"), crs = 4326)

trail <- rbind(trail_start, trail_end)

School$dist_to_park <- st_distance(School, Park) %>% apply(1, min) %>% as.numeric()
School$dist_to_hospital <- st_distance(School, Hospital) %>% apply(1, min) %>% as.numeric()
School$dist_to_clinc <- st_distance(School, clinic) %>% apply(1, min) %>% as.numeric()
School$dist_to_sports <- st_distance(School, Sports_Center) %>% apply(1, min) %>% as.numeric()
School$dist_to_health <- st_distance(School, health_food) %>% apply(1, min) %>% as.numeric()
School$dist_to_fastfood <- st_distance(School, fastfood) %>% apply(1, min) %>% as.numeric()
School$dist_to_river <- st_distance(School, river) %>% apply(1, min) %>% as.numeric()
School$dist_to_trail <- st_distance(School, trail) %>% apply(1, min) %>% as.numeric()

# 各區健康生活機能指數構建（將距離轉換為 0~1 的標準化「分數」）
normalize_score <- function(x) {
  1 - (x - min(x)) / (max(x) - min(x))
}

# 速食店距離標準化函數，距離越近分數越低（健康機能負面指標，扣分）
normalize_score_fastfood <- function(x) {
  (x - min(x)) / (max(x) - min(x))  # 距離越大分數越高（越健康）
  }

# 運動類指標標準化（距離越近越好）
School$score_park <- normalize_score(School$dist_to_park)
School$score_trail <- normalize_score(School$dist_to_trail)
School$score_sports <- normalize_score(School$dist_to_sports)
School$score_river <- normalize_score(School$dist_to_river)

# 飲食類指標標準化（健康飲食距離越近越好，速食距離越遠越好）
School$score_health <- normalize_score(School$dist_to_health)
School$score_fastfood <- normalize_score_fastfood(School$dist_to_fastfood)

# 醫療類指標標準化（距離越近越好）
School$score_hospital <- normalize_score(School$dist_to_hospital)
School$score_clinic <- normalize_score(School$dist_to_clinc)

# 合成三大類指數（簡單平均）
School$sport_index <- rowMeans(st_drop_geometry(School)[, c("score_park", "score_trail", "score_sports", "score_river")], na.rm = TRUE)
School$diet_index <- rowMeans(st_drop_geometry(School)[, c("score_health", "score_fastfood")], na.rm = TRUE)
School$medical_index <- rowMeans(st_drop_geometry(School)[, c("score_hospital", "score_clinic")], na.rm = TRUE)

# 計算整體健康生活機能指數（三項平均）
School$overall_health_index <- rowMeans(st_drop_geometry(School)[, c("sport_index", "diet_index", "medical_index")], na.rm = TRUE)

library(tmap)

tmap_mode("view")

map_sport <- tm_shape(School) +
  tm_dots(col = "sport_index", palette = "Blues", size = 0.5, title = "運動指數") +
  tm_layout(title = "台北市國小周邊運動生活機能指數", title.size = 1.2, title.position = c("center", "top"))

map_diet <- tm_shape(School) +
  tm_dots(col = "diet_index", palette = "Oranges", size = 0.5, title = "飲食指數") +
  tm_layout(title = "台北市國小周邊飲食生活機能指數", title.size = 1.2, title.position = c("center", "top"))

map_medical <- tm_shape(School) +
  tm_dots(col = "medical_index", palette = "Greens", size = 0.5, title = "醫療指數") +
  tm_layout(title = "台北市國小周邊醫療生活機能指數", title.size = 1.2, title.position = c("center", "top"))

map_overall <- tm_shape(School) +
  tm_dots(col = "overall_health_index", palette = "RdYlGn", size = 0.5, title = "整體健康指數") +
  tm_layout(title = "台北市國小周邊整體健康生活機能指數", title.size = 1.2, title.position = c("center", "top"))

tmap_arrange(map_sport, map_diet, map_medical, map_overall, ncol = 2)

NA

School %>%
  slice_min(order_by = sport_index, n = 10) %>%
  arrange(desc(sport_index)) %>%
  mutate(
    rank_factor = factor(Name, levels = Name),
    fill_color = gradient_n_pal(c("#0066CC", "#CCE5FF"))(rescale(sport_index))  # 深→淺
  ) %>%
  ggplot(aes(x = rank_factor, y = sport_index, fill = fill_color)) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  scale_fill_identity() +
  theme_minimal() +
  labs(
    x = "學校名稱",
    y = "指數分數",
    title = "臺北市國小周邊運動生活指數倒數 10 名"
  ) +
  theme(
    text = element_text(family = "PingFangSC-SemiBold"),
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 12, face = "bold")
  )

School %>%
  slice_min(order_by = diet_index, n = 10) %>%
  arrange(desc(diet_index)) %>%
  mutate(
    rank_factor = factor(Name, levels = Name),
    fill_color = gradient_n_pal(c("#FF6600", "#FFE5CC"))(rescale(diet_index))  # 深→淺
  ) %>%
  ggplot(aes(x = rank_factor, y = diet_index, fill = fill_color)) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  scale_fill_identity() +
  theme_minimal() +
  labs(
    x = "學校名稱",
    y = "指數分數",
    title = "臺北市國小周邊飲食生活指數倒數 10 名"
  ) +
  theme(
    text = element_text(family = "PingFangSC-SemiBold"),
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 12, face = "bold")
  )

School %>%
  slice_min(order_by = medical_index, n = 10) %>%
  arrange(desc(medical_index)) %>%
  mutate(
    rank_factor = factor(Name, levels = Name),
    fill_color = gradient_n_pal(c("#009933", "#CCFFCC"))(rescale(medical_index))  # 深→淺
  ) %>%
  ggplot(aes(x = rank_factor, y = medical_index, fill = fill_color)) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  scale_fill_identity() +
  theme_minimal() +
  labs(
    x = "學校名稱",
    y = "指數分數",
    title = "臺北市國小周邊醫療生活指數倒數 10 名"
  ) +
  theme(
    text = element_text(family = "PingFangSC-SemiBold"),
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 12, face = "bold"))

School %>%
  slice_min(order_by = overall_health_index, n = 10) %>%
  arrange(desc(overall_health_index)) %>%
  mutate(
    rank_factor = factor(Name, levels = Name),
    fill_color = gradient_n_pal(c("#D73027", "#FEE08B"))(rescale(overall_health_index))  # 深→淺
  ) %>%
  ggplot(aes(x = rank_factor, y = overall_health_index, fill = fill_color)) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  scale_fill_identity() +
  theme_minimal() +
  labs(
    x = "學校名稱",
    y = "指數分數",
    title = "臺北市國小周邊整體健康生活指數倒數 10 名"
  ) +
  theme(
    text = element_text(family = "PingFangSC-SemiBold"),
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 12, face = "bold"))

School_district <- st_join(School, taipei, left = TRUE)
# 移除空間資訊，並依行政區分組計算平均指數
district_summary <- School_district %>%
  st_set_geometry(NULL) %>%
  group_by(TOWNNAME) %>%
  summarise(
    sport_index_avg = mean(sport_index, na.rm = TRUE),
    diet_index_avg = mean(diet_index, na.rm = TRUE),
    medical_index_avg = mean(medical_index, na.rm = TRUE),
    overall_health_index_avg = mean(overall_health_index, na.rm = TRUE)
  )

# 將計算好的平均值合併回 taipei 空間資料
taipei <- taipei %>%
  left_join(district_summary, by = "TOWNNAME")

# 繪圖
ggplot(taipei) +
  geom_sf(aes(fill = overall_health_index_avg), color = "white") +
  scale_fill_gradientn(
  colours = c("#D73027", "#FEE08B", "#1A9850"),
  na.value = "grey90") +
  labs(fill = "整體健康指數平均", title = "台北市各行政區之整體健康指數") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), text = element_text(family = "PingFangSC-SemiBold"))

# 繪圖
ggplot(taipei) +
  geom_sf(aes(fill = sport_index_avg), color = "white") +
  scale_fill_gradientn(
  colours = c("#D73027", "#FEE08B", "#1A9850"),
  na.value = "grey90") +
  labs(fill = "運動健康指數平均", title = "台北市各行政區之運動健康指數") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), text = element_text(family = "PingFangSC-SemiBold"))

# 建立鄰接關係
nb <- poly2nb(taipei)

# 轉換為空間權重矩陣（row-standardized）
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# 計算 Local Moran's I
local_moran <- localmoran(taipei$overall_health_index_avg, lw, zero.policy = TRUE)

# 加入結果欄位
taipei$local_I <- local_moran[, "Ii"]
taipei$local_I_p <- local_moran[, "Pr(z != E(Ii))"]  # 或用自算的 p 值

# 畫圖：Local Moran's I 值分布
ggplot(taipei) +
  geom_sf(aes(fill = local_I), color = "white") +
  scale_fill_gradientn(
    colours = c("#D73027", "#FEE08B", "#1A9850"),
    na.value = "grey90"
  ) +
  labs(
    fill = "Local Moran's I",
    title = "台北市行政區之 Local Moran's I（整體健康指數）"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5),
    text = element_text(family = "PingFangSC-SemiBold")
  )