Spatial Analysis Final

Reading Data

rm(list = ls())

# The Taipei data provided in class
tpe_vill_sf = st_read("C:/113-2_Spring/Spatial_Analysis/coursefiles/Taipei_Vill/Taipei_Vill.shp", options = "ENCODING=big5")

options:        ENCODING=big5 
Reading layer `Taipei_Vill' from data source 
  `C:\113-2_Spring\Spatial_Analysis\coursefiles\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

# Population data
tpe_pop_df = read.csv("C:/113-2_Spring/Spatial_Analysis/final_project/0608_Plotting/taipei_vill_pop.csv", fileEncoding = "big5")

# All schools in Taiwan
schools <- read.csv("C:/113-2_Spring/Spatial_Analysis/final_project/0608_Plotting/school_112_csv.csv",
                     fileEncoding = "BIG5")

## Select those whose coordinate is not missing
schools_clean <- schools %>%
  filter(!is.na(`X.坐標`) & !is.na(`Y.坐標`))

## Turning into sf
schools_sf <- st_as_sf(schools_clean,
                       coords = c("X.坐標", "Y.坐標"),
                       crs = 3826)  # EPSG:3826 是 TWD97 二度分帶

# Sportfield Data
sportfield_sf = st_read("C:/113-2_Spring/Spatial_Analysis/final_project/0608_Plotting/sportfield_taipei/sportfield_taipei.shp")

Reading layer `sportfield_taipei' from data source 
  `C:\113-2_Spring\Spatial_Analysis\final_project\0608_Plotting\sportfield_taipei\sportfield_taipei.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1359 features and 51 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 121.4618 ymin: 24.97106 xmax: 121.6198 ymax: 25.1691
Geodetic CRS:  WGS 84

# Gaming places data
gaming_sf = st_read("C:/113-2_Spring/Spatial_Analysis/final_project/0608_Plotting/Gaming_center/Gaming_center.shp")

Reading layer `Gaming_center' from data source 
  `C:\113-2_Spring\Spatial_Analysis\final_project\0608_Plotting\Gaming_center\Gaming_center.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 105 features and 43 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 121.4959 ymin: 24.98938 xmax: 121.6208 ymax: 25.1241
Geodetic CRS:  WGS 84

# Hospital data
hospital_sf = st_read("C:/113-2_Spring/Spatial_Analysis/final_project/0608_Plotting/Hospital_taipei/Hospital_taipei.shp")

Reading layer `Hospital_taipei' from data source 
  `C:\113-2_Spring\Spatial_Analysis\final_project\0608_Plotting\Hospital_taipei\Hospital_taipei.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 876 features and 117 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 121.493 ymin: 24.98035 xmax: 121.6211 ymax: 25.1385
Geodetic CRS:  WGS 84

# fastfood data
fastfood_sf = st_read("C:/113-2_Spring/Spatial_Analysis/final_project/0608_Plotting/Tpe_Fastfood/Tpe_Fastfood.shp")

Reading layer `Tpe_Fastfood' from data source 
  `C:\113-2_Spring\Spatial_Analysis\final_project\0608_Plotting\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

# park data
park_sf = st_read("C:/113-2_Spring/Spatial_Analysis/final_project/0608_Plotting/Park_centroids/Park_centroids.shp")

Reading layer `Park_centroids' from data source 
  `C:\113-2_Spring\Spatial_Analysis\final_project\0608_Plotting\Park_centroids\Park_centroids.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 3115 features and 169 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 121.4622 ymin: 24.96698 xmax: 121.6249 ymax: 25.18077
Geodetic CRS:  WGS 84

# Selecting population data
tpe_pop_df = tpe_pop_df[tpe_pop_df$性別 == "計" & tpe_pop_df$區域代碼 > 100000000 & tpe_pop_df$年份 == 113 & tpe_pop_df$月份 == 2,  ]
tpe_vill_sf$區域代碼 = as.numeric(tpe_vill_sf$TOWN_ID) * 1000 + as.numeric(tpe_vill_sf$VILLAGE_ID)

Calculation

# Calculations on population 
tpe_pop_df$age_1_6 = rowSums(tpe_pop_df[8:13])
tpe_pop_df$age_7_12 = rowSums(tpe_pop_df[14:19])
tpe_pop_df$adult = rowSums(tpe_pop_df[26:71])

# primary_ratio and secondary ratio
tpe_pop_df$primary_ratio = tpe_pop_df$age_1_6 / tpe_pop_df$adult
tpe_pop_df$secondary_ratio = tpe_pop_df$age_7_12 / tpe_pop_df$adult

tpe_vill_sf = merge(tpe_vill_sf, tpe_pop_df, by = "區域代碼")

Selecting Schools

# 選取臺北的國中
schools_sf = schools_sf %>%
  filter(學校級別 %in% c("國民中學", "附設國民中學")) %>%
  filter(縣市別 == "臺北市")

# 選取「北市公立五虎 + 景美」
five_tigers_sf = schools_sf %>%
  filter(學校名稱 %in% c("市立中正國中", "市立龍門國中", "市立敦化國中", "市立介壽國中", "市立金華國中", "市立景美國中"))

# 選取其他建北率高的學校
other_top15_sf = schools_sf %>%
  filter(學校名稱 %in% c("國立師大附中附設國中", "私立靜心高中附設國中", "私立延平中學附設國中", "市立大安國中", "市立麗山國中","市立南門國中", "市立天母國中", "市立興雅國中", "市立石牌國中"))

# 建立簡化名稱欄位
five_tigers_sf$school_short <- recode(five_tigers_sf$學校名稱,
                                      "市立中正國中" = "中正",
                                      "市立龍門國中" = "龍門",
                                      "市立敦化國中" = "敦化",
                                      "市立介壽國中" = "介壽",
                                      "市立金華國中" = "金華",
                                      "市立景美國中" = "景美")

LISA of primary ratio (1~6)

nb = poly2nb(tpe_vill_sf, queen = FALSE) 
listw = nb2listw(nb, style = "W", zero.policy = TRUE)
density = as.vector(tpe_vill_sf$primary_ratio)
lisa = localmoran(density, listw = listw, zero.policy = TRUE)

tpe_vill_sf$local_I_pri = lisa[, 1]
tpe_vill_sf$z_score_pri = lisa[, 4]
tpe_vill_sf$p_value_pri = lisa[, 5]
quadrant = attr(lisa, "quadr")$mean
quadrant = factor(quadrant, levels = c(levels(quadrant), "NoSig"))
signif_level_05 = 0.05
quadrant[lisa[, 5] > signif_level_05] = "NoSig"
tpe_vill_sf$LISA_type_pri = quadrant

lisa_colors <- c(
  'High-High' = 'red',
  'Low-Low' = 'blue',
  'High-Low' = 'pink',
  'Low-High' = 'skyblue2',
  'NoSig' = 'grey'
)

tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_pri",
              palette = lisa_colors,
              title = "LISA Cluster Type (age 1~6 / adult)") +
  tm_layout(title = "LISA Map for Age 1~6 / Adult", legend.outside = TRUE)

LISA of secondary ratio (7~12)

nb = poly2nb(tpe_vill_sf, queen = FALSE) 
listw = nb2listw(nb, style = "W", zero.policy = TRUE)
density = as.vector(tpe_vill_sf$secondary_ratio)
lisa = localmoran(density, listw = listw, zero.policy = TRUE)

tpe_vill_sf$local_I_sec = lisa[, 1]
tpe_vill_sf$z_score_sec = lisa[, 4]
tpe_vill_sf$p_value_sec = lisa[, 5]
quadrant = attr(lisa, "quadr")$mean
quadrant = factor(quadrant, levels = c(levels(quadrant), "NoSig"))
signif_level = 0.05
quadrant[lisa[, 5] > signif_level] = "NoSig"
tpe_vill_sf$LISA_type_sec = quadrant

tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_sec",
              palette = lisa_colors,
              title = "LISA Cluster Type (age 7~12 / adult)") +
  tm_layout(title = "LISA Map for Age 7~12 / Adult", legend.outside = TRUE)

density_sec = as.vector(tpe_vill_sf$secondary_ratio)
lisa_sec = localmoran(density_sec, listw = listw, zero.policy = TRUE)

tpe_vill_sf$local_I_sec = lisa_sec[, 1]
tpe_vill_sf$z_score_sec = lisa_sec[, 4]
tpe_vill_sf$p_value_sec = lisa_sec[, 5]

quadrant_sec = attr(lisa_sec, "quadr")$mean
quadrant_sec = factor(quadrant_sec, levels = c(levels(quadrant_sec), "NoSig"))

# 使用 FDR 校正
quadrant_sec[p.adjust(lisa_sec[, 5], method = "fdr") > signif_level] = "NoSig"
tpe_vill_sf$LISA_type_sec = quadrant_sec

# 繪圖（含學校點位）
tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_sec",
              palette = lisa_colors,
              title = "LISA Cluster Type (age 7~12 / adult)") +
  tm_layout(title = "LISA Map for Age 7~12 / Adult", 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(title = )`

Displaying the five tigers on the FDRd map

names(grDevices::windowsFonts())  # Windows 專用

[1] "serif"        "sans"         "mono"         "Noto Sans TC"

# 手動為每個學校加上偏移值（你可以依照實際需求微調）
five_tigers_sf <- five_tigers_sf %>%
  mutate(
    xmod = c(0.3, -0.5, 0, 1, -1, 1),   # X 偏移
    ymod = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8) # Y 偏移
  )

tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_sec",
              palette = lisa_colors,
              title = "LISA Cluster Type (age 7~12 / adult)") +
  
  # 其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +
  
  # 公立五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +
  
  # 主文字
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod",
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +
  
  tm_layout(title = "LISA Map for Age 7~12 / Adult",
            legend.outside = TRUE)

Environment

Gaming Center

# Setting sig_level
signif_gaming = 0.01

gaming_sf <- st_transform(gaming_sf, st_crs(tpe_vill_sf))

joined_gaming <- st_join(gaming_sf, tpe_vill_sf, join = st_within)

gaming_counts <- joined_gaming %>%
  group_by(TOWN_ID) %>%
  summarise(gaming_count = n())

gaming_counts_df <- st_drop_geometry(gaming_counts)

tpe_vill_sf$gaming_count <- gaming_counts_df$gaming_count[match(tpe_vill_sf$TOWN_ID, gaming_counts_df$TOWN_ID)]
tpe_vill_sf$gaming_count[is.na(tpe_vill_sf$gaming_count)] <- 0

tpe_vill_sf$AREA <- set_units(st_area(tpe_vill_sf), "km^2")
tpe_vill_sf$gaming_density <- tpe_vill_sf$gaming_count / tpe_vill_sf$AREA

coords <- st_coordinates(st_centroid(tpe_vill_sf))

nb <- dnearneigh(coords, d1 = 0, d2 = 10000)

dist_list <- nbdists(nb, coords)

inv_dist_list <- lapply(dist_list, function(x) {
  y <- 1 / (x + 0.0001)
  return(y)
})

listw <- nb2listw(nb, glist = inv_dist_list, style = "W", zero.policy = TRUE)

density <- as.vector(tpe_vill_sf$gaming_density)

global_moran <- moran.test(density, listw = listw, zero.policy = TRUE)
print(global_moran)

    Moran I test under randomisation

data:  density  
weights: listw    

Moran I statistic standard deviate = 46.138, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     1.775105e-01     -2.197802e-03      1.517092e-05 

lisa <- localmoran(density, listw = listw, zero.policy = TRUE)

tpe_vill_sf$local_I <- lisa[, 1]
tpe_vill_sf$z_score <- lisa[, 4]
tpe_vill_sf$p_value <- lisa[, 5]

quadrant <- attr(lisa, "quadr")$mean
quadrant <- factor(quadrant, levels = c(levels(quadrant), "NoSig"))

tpe_vill_sf$pvalue_FDR <- p.adjust(lisa[, 5], method = "fdr")

tpe_vill_sf$LISA_type_Gaming <- "NoSig"
tpe_vill_sf$LISA_type_Gaming[tpe_vill_sf$pvalue_FDR < signif_gaming & tpe_vill_sf$local_I > 0 & quadrant == "High-High"] <- "High-High"
tpe_vill_sf$LISA_type_Gaming[tpe_vill_sf$pvalue_FDR < signif_gaming & tpe_vill_sf$local_I > 0 & quadrant == "Low-Low"] <- "Low-Low"
tpe_vill_sf$LISA_type_Gaming[tpe_vill_sf$pvalue_FDR < signif_gaming & tpe_vill_sf$local_I < 0 & quadrant == "High-Low"] <- "High-Low"
tpe_vill_sf$LISA_type_Gaming[tpe_vill_sf$pvalue_FDR < signif_gaming & tpe_vill_sf$local_I < 0 & quadrant == "Low-High"] <- "Low-High"

tpe_vill_sf$LISA_type_Gaming <- factor(tpe_vill_sf$LISA_type_Gaming,
  levels = c("High-High", "Low-Low", "High-Low", "Low-High", "NoSig"))

lisa_colors <- c(
  'High-High' = 'red',
  'Low-Low' = 'blue',
  'High-Low' = 'pink',
  'Low-High' = 'skyblue2',
  'NoSig' = 'grey'
)

tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_Gaming",
              palette = lisa_colors,
              title = "LISA Cluster Type (Sport Facilities)") +
  
  # 其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +
  
  # 五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +
  
  # 主文字
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod",
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +
  
  tm_layout(title = "LISA Map (Inverse Distance) - Gaming Centers",
            legend.outside = TRUE)

fastfood

# 設定信心水準
signif_fastfood <- 0.05

# 投影一致
fastfood_sf <- st_transform(fastfood_sf, st_crs(tpe_vill_sf))

# 空間連接：找出每個村里有幾個 fast food 點
joined_fastfood <- st_join(fastfood_sf, tpe_vill_sf, join = st_within)

fastfood_counts <- joined_fastfood %>%
  group_by(TOWN_ID) %>%
  summarise(fastfood_count = n())

fastfood_counts_df <- st_drop_geometry(fastfood_counts)

tpe_vill_sf$fastfood_count <- fastfood_counts_df$fastfood_count[match(tpe_vill_sf$TOWN_ID, fastfood_counts_df$TOWN_ID)]
tpe_vill_sf$fastfood_count[is.na(tpe_vill_sf$fastfood_count)] <- 0

# 面積與密度
tpe_vill_sf$AREA <- set_units(st_area(tpe_vill_sf), "km^2")
tpe_vill_sf$fastfood_density <- tpe_vill_sf$fastfood_count / tpe_vill_sf$AREA

# 鄰近矩陣設定
coords <- st_coordinates(st_centroid(tpe_vill_sf))

警告： st_centroid assumes attributes are constant over geometries

nb <- dnearneigh(coords, d1 = 0, d2 = 10000)

dist_list <- nbdists(nb, coords)
inv_dist_list <- lapply(dist_list, function(x) 1 / (x + 0.0001))

listw <- nb2listw(nb, glist = inv_dist_list, style = "W", zero.policy = TRUE)

# 執行 LISA
density <- as.vector(tpe_vill_sf$fastfood_density)

global_moran <- moran.test(density, listw = listw, zero.policy = TRUE)
print(global_moran)

    Moran I test under randomisation

data:  density  
weights: listw    

Moran I statistic standard deviate = 39.132, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     1.529806e-01     -2.197802e-03      1.572547e-05 

lisa <- localmoran(density, listw = listw, zero.policy = TRUE)

tpe_vill_sf$local_I <- lisa[, 1]
tpe_vill_sf$z_score <- lisa[, 4]
tpe_vill_sf$p_value <- lisa[, 5]

quadrant <- attr(lisa, "quadr")$mean
quadrant <- factor(quadrant, levels = c(levels(quadrant), "NoSig"))

tpe_vill_sf$pvalue_FDR <- p.adjust(lisa[, 5], method = "fdr")

tpe_vill_sf$LISA_type_Fastfood <- "NoSig"
tpe_vill_sf$LISA_type_Fastfood[tpe_vill_sf$pvalue_FDR < signif_fastfood & tpe_vill_sf$local_I > 0 & quadrant == "High-High"] <- "High-High"
tpe_vill_sf$LISA_type_Fastfood[tpe_vill_sf$pvalue_FDR < signif_fastfood & tpe_vill_sf$local_I > 0 & quadrant == "Low-Low"] <- "Low-Low"
tpe_vill_sf$LISA_type_Fastfood[tpe_vill_sf$pvalue_FDR < signif_fastfood & tpe_vill_sf$local_I < 0 & quadrant == "High-Low"] <- "High-Low"
tpe_vill_sf$LISA_type_Fastfood[tpe_vill_sf$pvalue_FDR < signif_fastfood & tpe_vill_sf$local_I < 0 & quadrant == "Low-High"] <- "Low-High"

tpe_vill_sf$LISA_type_Fastfood <- factor(tpe_vill_sf$LISA_type_Fastfood,
                                         levels = c("High-High", "Low-Low", "High-Low", "Low-High", "NoSig"))

# 配色
lisa_colors <- c(
  'High-High' = 'red',
  'Low-Low' = 'blue',
  'High-Low' = 'pink',
  'Low-High' = 'skyblue2',
  'NoSig' = 'grey'
)

# 畫出地圖
tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_Fastfood",
              palette = lisa_colors,
              title = "LISA Cluster Type (Fast Food)") +
  
  # 其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +
  
  # 五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +
  
  # 標上學校名稱
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod",
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +
  
  tm_layout(title = "LISA Map (Inverse Distance) - Fast Food",
            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(title = )`

intersecting

# 建立交集標記欄位
tpe_vill_sf$HH_Intersect <- ifelse(
  tpe_vill_sf$LISA_type_Fastfood == "High-High" & 
    tpe_vill_sf$LISA_type_Gaming == "High-High", 
  "High-High Both", 
  "Other"
)

# 設定顏色
intersect_colors <- c("High-High Both" = "purple", "Other" = "lightgrey")

# 畫出地圖
tm_shape(tpe_vill_sf) +
  tm_polygons("HH_Intersect", 
              palette = intersect_colors, 
              title = "High-High Intersection") +
  
  # 其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +
  
  # 五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +
  
  # 標上學校名稱
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod",
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +

  tm_layout(title = "Intersection of High-High Clusters: Fast Food & Gaming Centers",
            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(title = )`

sportfield

sportfield_sf <- st_transform(sportfield_sf, st_crs(tpe_vill_sf))

joined_sportfield <- st_join(sportfield_sf, tpe_vill_sf, join = st_within)

sportfield_counts <- joined_sportfield %>%
  group_by(TOWN_ID) %>%
  summarise(sportfield_count = n())

sportfield_counts_df <- st_drop_geometry(sportfield_counts)

tpe_vill_sf$sportfield_count <- sportfield_counts_df$sportfield_count[match(tpe_vill_sf$TOWN_ID, sportfield_counts_df$TOWN_ID)]
tpe_vill_sf$sportfield_count[is.na(tpe_vill_sf$sportfield_count)] <- 0

tpe_vill_sf$AREA <- set_units(st_area(tpe_vill_sf), "km^2")
tpe_vill_sf$sportfield_density <- tpe_vill_sf$sportfield_count / tpe_vill_sf$AREA

coords <- st_coordinates(st_centroid(tpe_vill_sf))

# 10km
nb <- dnearneigh(coords, d1 = 0, d2 = 10000)

dist_list <- nbdists(nb, coords)

inv_dist_list <- lapply(dist_list, function(x) {
  y <- 1 / (x + 0.0001)
  return(y)
})

listw <- nb2listw(nb, glist = inv_dist_list, style = "W", zero.policy = TRUE)

density <- as.vector(tpe_vill_sf$sportfield_density)

global_moran <- moran.test(density, listw = listw, zero.policy = TRUE)
print(global_moran)

    Moran I test under randomisation

data:  density  
weights: listw    

Moran I statistic standard deviate = 22.131, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     8.543080e-02     -2.197802e-03      1.567804e-05 

lisa <- localmoran(density, listw = listw, zero.policy = TRUE)

tpe_vill_sf$local_I <- lisa[, 1]
tpe_vill_sf$z_score <- lisa[, 4]
tpe_vill_sf$p_value <- lisa[, 5]

quadrant <- attr(lisa, "quadr")$mean
quadrant <- factor(quadrant, levels = c(levels(quadrant), "NoSig"))

tpe_vill_sf$pvalue_FDR <- p.adjust(lisa[, 5], method = "fdr")

# Setting sig_level
signif_sportfield = 0.01

tpe_vill_sf$LISA_type_sportfield <- "NoSig"
tpe_vill_sf$LISA_type_sportfield[tpe_vill_sf$pvalue_FDR < signif_sportfield & tpe_vill_sf$local_I > 0 & quadrant == "High-High"] <- "High-High"
tpe_vill_sf$LISA_type_sportfield[tpe_vill_sf$pvalue_FDR < signif_sportfield & tpe_vill_sf$local_I > 0 & quadrant == "Low-Low"] <- "Low-Low"
tpe_vill_sf$LISA_type_sportfield[tpe_vill_sf$pvalue_FDR < signif_sportfield & tpe_vill_sf$local_I < 0 & quadrant == "High-Low"] <- "High-Low"
tpe_vill_sf$LISA_type_sportfield[tpe_vill_sf$pvalue_FDR < signif_sportfield & tpe_vill_sf$local_I < 0 & quadrant == "Low-High"] <- "Low-High"

tpe_vill_sf$LISA_type_sportfield <- factor(tpe_vill_sf$LISA_type_sportfield,
  levels = c("High-High", "Low-Low", "High-Low", "Low-High", "NoSig"))

lisa_colors <- c(
  'High-High' = 'red',
  'Low-Low' = 'blue',
  'High-Low' = 'pink',
  'Low-High' = 'skyblue2',
  'NoSig' = 'grey'
)

tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_sportfield",
              palette = lisa_colors,
              title = "LISA Cluster Type (Sport Facilities)") +
  
  # 其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +
  
  # 五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +
  
  # 主文字
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod",
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +
  
  tm_layout(title = "LISA Map (Inverse Distance) - Sport Facilities",
            legend.outside = TRUE)

Hospital

# 設定信心水準
signif_hospital <- 0.01

# 確保 hospital 資料與行政區使用同樣 CRS
hospital_sf <- st_transform(hospital_sf, st_crs(tpe_vill_sf))

# 加入行政區
joined_hospital <- st_join(hospital_sf, tpe_vill_sf, join = st_within)

# 每個行政區的醫療點數量
hospital_counts <- joined_hospital %>%
  group_by(TOWN_ID) %>%
  summarise(hospital_count = n())

# 移除 geometry
hospital_counts_df <- st_drop_geometry(hospital_counts)

# 對應進原本的行政區圖層
tpe_vill_sf$hospital_count <- hospital_counts_df$hospital_count[match(tpe_vill_sf$TOWN_ID, hospital_counts_df$TOWN_ID)]
tpe_vill_sf$hospital_count[is.na(tpe_vill_sf$hospital_count)] <- 0

# 計算密度（設施數 / 區域面積）
tpe_vill_sf$AREA <- set_units(st_area(tpe_vill_sf), "km^2")
tpe_vill_sf$hospital_density <- tpe_vill_sf$hospital_count / tpe_vill_sf$AREA

# 計算鄰近矩陣
coords <- st_coordinates(st_centroid(tpe_vill_sf))
nb <- dnearneigh(coords, d1 = 0, d2 = 10000)

dist_list <- nbdists(nb, coords)

inv_dist_list <- lapply(dist_list, function(x) {
  1 / (x + 0.0001)
})

listw <- nb2listw(nb, glist = inv_dist_list, style = "W", zero.policy = TRUE)

# LISA 分析
density <- as.vector(tpe_vill_sf$hospital_density)

global_moran <- moran.test(density, listw = listw, zero.policy = TRUE)
print(global_moran)

    Moran I test under randomisation

data:  density  
weights: listw    

Moran I statistic standard deviate = 40.53, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     1.585062e-01     -2.197802e-03      1.572152e-05 

lisa <- localmoran(density, listw = listw, zero.policy = TRUE)

tpe_vill_sf$local_I <- lisa[, 1]
tpe_vill_sf$z_score <- lisa[, 4]
tpe_vill_sf$p_value <- lisa[, 5]

quadrant <- attr(lisa, "quadr")$mean
quadrant <- factor(quadrant, levels = c(levels(quadrant), "NoSig"))

tpe_vill_sf$pvalue_FDR <- p.adjust(lisa[, 5], method = "fdr")

# 指定類型
tpe_vill_sf$LISA_type_Hospital <- "NoSig"
tpe_vill_sf$LISA_type_Hospital[tpe_vill_sf$pvalue_FDR < signif_hospital & tpe_vill_sf$local_I > 0 & quadrant == "High-High"] <- "High-High"
tpe_vill_sf$LISA_type_Hospital[tpe_vill_sf$pvalue_FDR < signif_hospital & tpe_vill_sf$local_I > 0 & quadrant == "Low-Low"] <- "Low-Low"
tpe_vill_sf$LISA_type_Hospital[tpe_vill_sf$pvalue_FDR < signif_hospital & tpe_vill_sf$local_I < 0 & quadrant == "High-Low"] <- "High-Low"
tpe_vill_sf$LISA_type_Hospital[tpe_vill_sf$pvalue_FDR < signif_hospital & tpe_vill_sf$local_I < 0 & quadrant == "Low-High"] <- "Low-High"

tpe_vill_sf$LISA_type_Hospital <- factor(tpe_vill_sf$LISA_type_Hospital,
                                         levels = c("High-High", "Low-Low", "High-Low", "Low-High", "NoSig"))

# 配色
lisa_colors <- c(
  'High-High' = 'red',
  'Low-Low' = 'blue',
  'High-Low' = 'pink',
  'Low-High' = 'skyblue2',
  'NoSig' = 'grey'
)

# 畫圖
tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_Hospital",
              palette = lisa_colors,
              title = "LISA Cluster Type (Hospitals)") +
  
  # 其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +
  
  # 五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +
  
  # 主文字
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod",
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +
  
  tm_layout(title = "LISA Map (Inverse Distance) - Hospitals",
            legend.outside = TRUE)

Park

# 1. 設定信心水準
signif_park <- 0.01

# 2. 座標轉換
park_sf <- st_transform(park_sf, st_crs(tpe_vill_sf))

# 3. Spatial join 與計數
joined_park <- st_join(park_sf, tpe_vill_sf, join = st_within)

park_counts <- joined_park %>%
  group_by(TOWN_ID) %>%
  summarise(park_count = n())

park_counts_df <- st_drop_geometry(park_counts)

# 4. 加入計數與密度欄位
tpe_vill_sf$park_count <- park_counts_df$park_count[match(tpe_vill_sf$TOWN_ID, park_counts_df$TOWN_ID)]
tpe_vill_sf$park_count[is.na(tpe_vill_sf$park_count)] <- 0

tpe_vill_sf$AREA <- set_units(st_area(tpe_vill_sf), "km^2")
tpe_vill_sf$park_density <- tpe_vill_sf$park_count / tpe_vill_sf$AREA

# 5. 計算鄰近矩陣與空間權重
coords <- st_coordinates(st_centroid(tpe_vill_sf))

nb <- dnearneigh(coords, d1 = 0, d2 = 10000)
dist_list <- nbdists(nb, coords)

inv_dist_list <- lapply(dist_list, function(x) {
  1 / (x + 0.0001)
})

listw <- nb2listw(nb, glist = inv_dist_list, style = "W", zero.policy = TRUE)

# 6. Global & Local Moran's I
density <- as.vector(tpe_vill_sf$park_density)

global_moran <- moran.test(density, listw = listw, zero.policy = TRUE)
print(global_moran)

    Moran I test under randomisation

data:  density  
weights: listw    

Moran I statistic standard deviate = 52.42, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     2.048055e-01     -2.197802e-03      1.559429e-05 

lisa <- localmoran(density, listw = listw, zero.policy = TRUE)

# 7. 加入 LISA 指標
tpe_vill_sf$local_I_park <- lisa[, 1]
tpe_vill_sf$z_score_park <- lisa[, 4]
tpe_vill_sf$p_value_park <- lisa[, 5]

quadrant_park <- attr(lisa, "quadr")$mean
quadrant_park <- factor(quadrant_park, levels = c(levels(quadrant_park), "NoSig"))

# FDR 校正
tpe_vill_sf$pvalue_FDR_park <- p.adjust(tpe_vill_sf$p_value_park, method = "fdr")

# 分群
tpe_vill_sf$LISA_type_Park <- "NoSig"
tpe_vill_sf$LISA_type_Park[tpe_vill_sf$pvalue_FDR_park < signif_park & tpe_vill_sf$local_I_park > 0 & quadrant_park == "High-High"] <- "High-High"
tpe_vill_sf$LISA_type_Park[tpe_vill_sf$pvalue_FDR_park < signif_park & tpe_vill_sf$local_I_park > 0 & quadrant_park == "Low-Low"] <- "Low-Low"
tpe_vill_sf$LISA_type_Park[tpe_vill_sf$pvalue_FDR_park < signif_park & tpe_vill_sf$local_I_park < 0 & quadrant_park == "High-Low"] <- "High-Low"
tpe_vill_sf$LISA_type_Park[tpe_vill_sf$pvalue_FDR_park < signif_park & tpe_vill_sf$local_I_park < 0 & quadrant_park == "Low-High"] <- "Low-High"

tpe_vill_sf$LISA_type_Park <- factor(tpe_vill_sf$LISA_type_Park,
                                     levels = c("High-High", "Low-Low", "High-Low", "Low-High", "NoSig"))

# 畫圖
tm_shape(tpe_vill_sf) +
  tm_polygons("LISA_type_Park",
              palette = lisa_colors,
              title = "LISA Cluster Type (Parks)") +
  
  # 其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +
  
  # 五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +
  
  # 主文字
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod",
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +
  
  tm_layout(title = "LISA Map (Inverse Distance) - Parks",
            legend.outside = TRUE)

intersection

# 建立 High-High 與 Low-Low 的交集欄位
tpe_vill_sf$LISA_HHH_intersection <- with(tpe_vill_sf,
  LISA_type_Park == "High-High" &
  LISA_type_sportfield == "High-High" &
  LISA_type_Hospital == "High-High"
)

tpe_vill_sf$LISA_LLL_intersection <- with(tpe_vill_sf,
  LISA_type_Park == "Low-Low" &
  LISA_type_sportfield == "Low-Low" &
  LISA_type_Hospital == "Low-Low"
)

# 建立一個新的欄位標示交集類型
tpe_vill_sf$LISA_intersection_type <- "Other"
tpe_vill_sf$LISA_intersection_type[tpe_vill_sf$LISA_HHH_intersection] <- "All High-High"
tpe_vill_sf$LISA_intersection_type[tpe_vill_sf$LISA_LLL_intersection] <- "All Low-Low"
tpe_vill_sf$LISA_intersection_type <- factor(
  tpe_vill_sf$LISA_intersection_type,
  levels = c("Other", "All High-High", "All Low-Low")
)

# 畫圖
tm_shape(tpe_vill_sf) +
  # 畫出交集區域
  tm_polygons("LISA_intersection_type",
              palette = c("gray", "green", "darkred"),
              labels = c("Not all same type", "All High-High", "All Low-Low"),
              title = "LISA Intersection Type") +

  # 疊上其他建北率高的學校
  tm_shape(other_top15_sf) +
  tm_dots(size = 0.5, col = "gray", shape = 21, border.col = "white") +

  # 疊上五虎
  tm_shape(five_tigers_sf) +
  tm_dots(size = 0.5, col = "black", shape = 21, border.col = "white") +

  # 加上學校文字標籤
  tm_text("school_short", 
          size = 0.9, 
          xmod = "xmod", 
          ymod = "ymod", 
          col = "black", 
          fontfamily = "Noto Sans TC",
          fontface = "bold") +

  # 標題與圖例配置
  tm_layout(title = "Intersection of LISA Clusters (High-High and Low-Low)",
            legend.outside = TRUE)