---
disqus: ahb0222
GA : G-VF9ZT413CG
---
> [color=#40f1ef][name=LHB阿好伯][:earth_africa:](https://www.facebook.com/LHB0222/)
>[time=Tue, Dec 16, 2025 8:39 AM]

[TOC]

在環境監測專案中，我們經常需要回答一個問題

「某個地點周圍有多少監測站？」這支 R 腳本正是為了解決這個需求而設計。

核心功能： 透過 Haversine 公式計算球面距離
將測站依據與目標點的距離分為三個等級，並以 Leaflet 繪製互動式地圖。

三圈分類邏輯：

🔴 紅色標記：內圈範圍內（預設 1 公里）
🔵 藍色標記：內圈外、外圈內（預設 1-2 公里）
🟢 綠色標記：外圈範圍外（超過 2 公里）

彈性設計： 只需調整程式開頭的 radius_inner 和 radius_outer 兩個變數

即可自訂篩選半徑，適用於不同場景的分析需求。

地圖呈現特色：
===
支援 OpenStreetMap、CartoDB、衛星影像三種底圖切換
內圈紅色實線、外圈藍色虛線清楚標示範圍
點擊標記可查看測站名稱、位置及精確距離
圖例自動統計各類別數量

這個工具特別適合用於選址評估、監測網絡覆蓋分析，或是快速掌握特定區域的空品感測器分布狀況。


---

# 繪製測站地圖與半徑圓
# 使用 leaflet 建立互動式地圖

library(readr)
library(leaflet)
library(htmlwidgets)

# ============================================
# 可調整參數
# ============================================
radius_inner <- 1  # 內圈半徑（公里）- 紅色標記
radius_outer <- 2  # 外圈半徑（公里）- 藍色標記範圍
# 超過外圈半徑的為綠色標記

# 目標座標
target_lat <- 22.696152701092515
target_lon <- 120.2971322908144

# ============================================
# 讀取資料
# ============================================
df <- read_csv("20251212-0644.csv")

# Haversine 公式計算距離（公里）
haversine_distance <- function(lat1, lon1, lat2, lon2) {
  R <- 6371
  lat1_rad <- lat1 * pi / 180
  lat2_rad <- lat2 * pi / 180
  delta_lat <- (lat2 - lat1) * pi / 180
  delta_lon <- (lon2 - lon1) * pi / 180
  a <- sin(delta_lat / 2)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(delta_lon / 2)^2
  c <- 2 * atan2(sqrt(a), sqrt(1 - a))
  return(R * c)
}

# 計算每個測站與目標點的距離
df$distance_km <- mapply(
  haversine_distance,
  target_lat, target_lon,
  df$latitude, df$longitude
)

# ============================================
# 分類測站（三種類別）
# ============================================
df$category <- ifelse(
  df$distance_km <= radius_inner, 
  "inner",      # 紅色：內圈內
  ifelse(
    df$distance_km <= radius_outer, 
    "middle",   # 藍色：內圈外、外圈內
    "outer"     # 綠色：外圈外
  )
)

# 統計
stations_inner <- sum(df$category == "inner")
stations_middle <- sum(df$category == "middle")
stations_outer <- sum(df$category == "outer")

cat("============================================\n")
cat("搜尋中心: (", target_lat, ", ", target_lon, ")\n", sep = "")
cat("內圈半徑:", radius_inner, "公里\n")
cat("外圈半徑:", radius_outer, "公里\n")
cat("============================================\n")
cat("🔴 ", radius_inner, "公里內測站數量: ", stations_inner, "\n", sep = "")
cat("🔵 ", radius_inner, "-", radius_outer, "公里內測站數量: ", stations_middle, "\n", sep = "")
cat("🟢 ", radius_outer, "公里外測站數量: ", stations_outer, "\n", sep = "")
cat("============================================\n")

# ============================================
# 建立 leaflet 地圖
# ============================================
map <- leaflet() %>%
  # 底圖
  addTiles(group = "OpenStreetMap") %>%
  addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "衛星影像") %>%
  
  # 設定視圖中心為目標座標
  setView(lng = target_lon, lat = target_lat, zoom = 13) %>%
  
  # 繪製外圈半徑圓（藍色邊界）
  addCircles(
    lng = target_lon,
    lat = target_lat,
    radius = radius_outer * 1000,  # 轉換為公尺
    color = "#007BFF",
    weight = 2,
    fillColor = "#007BFF",
    fillOpacity = 0.05,
    dashArray = "5,5",
    popup = paste0(
      "<b>外圈範圍</b><br>",
      "半徑: ", radius_outer, " 公里<br>",
      "直徑: ", radius_outer * 2, " 公里<br>",
      "範圍內測站: ", stations_inner + stations_middle, " 個"
    ),
    group = "外圈範圍"
  ) %>%
  
  # 繪製內圈半徑圓（紅色邊界）
  addCircles(
    lng = target_lon,
    lat = target_lat,
    radius = radius_inner * 1000,  # 轉換為公尺
    color = "#DC3545",
    weight = 3,
    fillColor = "#DC3545",
    fillOpacity = 0.1,
    popup = paste0(
      "<b>內圈範圍</b><br>",
      "半徑: ", radius_inner, " 公里<br>",
      "直徑: ", radius_inner * 2, " 公里<br>",
      "範圍內測站: ", stations_inner, " 個"
    ),
    group = "內圈範圍"
  ) %>%
  
  # 標記目標中心點
  addMarkers(
    lng = target_lon,
    lat = target_lat,
    popup = paste0(
      "<b>🎯 搜尋中心點</b><br>",
      "緯度: ", round(target_lat, 6), "<br>",
      "經度: ", round(target_lon, 6), "<br>",
      "內圈: ", radius_inner, " 公里<br>",
      "外圈: ", radius_outer, " 公里"
    ),
    icon = makeIcon(
      iconUrl = "https://raw.githubusercontent.com/pointhi/leaflet-color-markers/master/img/marker-icon-red.png",
      iconWidth = 25, iconHeight = 41,
      iconAnchorX = 12, iconAnchorY = 41
    ),
    group = "中心點"
  ) %>%
  
  # 繪製內圈測站（紅色）
  addCircleMarkers(
    data = df[df$category == "inner", ],
    lng = ~longitude,
    lat = ~latitude,
    radius = 9,
    color = "#DC3545",
    fillColor = "#DC3545",
    fillOpacity = 0.9,
    weight = 2,
    popup = ~paste0(
      "<b>", thing_name, "</b><br>",
      "城市: ", city, "<br>",
      "鄉鎮: ", township, "<br>",
      "座標: (", round(latitude, 6), ", ", round(longitude, 6), ")<br>",
      "<span style='color: #DC3545; font-weight: bold;'>",
      "🔴 距離: ", round(distance_km, 3), " 公里 (", radius_inner, "km內)</span>"
    ),
    group = paste0(radius_inner, "km內測站")
  ) %>%
  
  # 繪製中圈測站（藍色）
  addCircleMarkers(
    data = df[df$category == "middle", ],
    lng = ~longitude,
    lat = ~latitude,
    radius = 7,
    color = "#007BFF",
    fillColor = "#007BFF",
    fillOpacity = 0.8,
    weight = 2,
    popup = ~paste0(
      "<b>", thing_name, "</b><br>",
      "城市: ", city, "<br>",
      "鄉鎮: ", township, "<br>",
      "座標: (", round(latitude, 6), ", ", round(longitude, 6), ")<br>",
      "<span style='color: #007BFF; font-weight: bold;'>",
      "🔵 距離: ", round(distance_km, 3), " 公里 (", radius_inner, "-", radius_outer, "km)</span>"
    ),
    group = paste0(radius_inner, "-", radius_outer, "km測站")
  ) %>%
  
  # 繪製外圈測站（綠色）
  addCircleMarkers(
    data = df[df$category == "outer", ],
    lng = ~longitude,
    lat = ~latitude,
    radius = 5,
    color = "#28A745",
    fillColor = "#28A745",
    fillOpacity = 0.6,
    weight = 1,
    popup = ~paste0(
      "<b>", thing_name, "</b><br>",
      "城市: ", city, "<br>",
      "鄉鎮: ", township, "<br>",
      "座標: (", round(latitude, 6), ", ", round(longitude, 6), ")<br>",
      "<span style='color: #28A745;'>",
      "🟢 距離: ", round(distance_km, 3), " 公里 (", radius_outer, "km外)</span>"
    ),
    group = paste0(radius_outer, "km外測站")
  ) %>%
  
  # 圖層控制
  addLayersControl(
    baseGroups = c("OpenStreetMap", "CartoDB", "衛星影像"),
    overlayGroups = c(
      "內圈範圍", 
      "外圈範圍", 
      "中心點", 
      paste0(radius_inner, "km內測站"),
      paste0(radius_inner, "-", radius_outer, "km測站"),
      paste0(radius_outer, "km外測站")
    ),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  
  # 圖例
  addLegend(
    position = "bottomright",
    colors = c("#DC3545", "#007BFF", "#28A745"),
    labels = c(
      paste0("🔴 ", radius_inner, "km內 (", stations_inner, ")"),
      paste0("🔵 ", radius_inner, "-", radius_outer, "km (", stations_middle, ")"),
      paste0("🟢 ", radius_outer, "km外 (", stations_outer, ")")
    ),
    title = paste0("測站分布 (共", nrow(df), "站)"),
    opacity = 0.9
  ) %>%
  
  # 比例尺
  addScaleBar(position = "bottomleft", options = scaleBarOptions(imperial = FALSE))

# 儲存為 HTML 檔案
saveWidget(map, "stations_map.html", selfcontained = TRUE)

cat("\n✅ 地圖已儲存為 stations_map.html\n")
cat("總測站數:", nrow(df), "\n")

# 多點位計算


# ============================================
# R語言實戰：用 Leaflet 打造空品測站距離篩選地圖
# 支援多點位清單，計算不重複站點數量
# ============================================

library(readr)
library(leaflet)
library(htmlwidgets)
library(dplyr)

# ============================================
# 可調整參數
# ============================================
radius_inner <- 1  # 內圈半徑（公里）- 紅色標記
radius_outer <- 2  # 外圈半徑（公里）- 藍色標記範圍

# ============================================
# 多點位清單設定
# 格式：data.frame(name, lat, lon)
# ============================================
target_points <- data.frame(
  name = c("點位A-鼓山區", "點位B-前鎮區", "點位C-左營區"),
  lat = c(22.696152701092515, 22.681366159432702, 22.690000),
  lon = c(120.2971322908144, 120.3790205926145, 120.295000)
)

# 點位顏色（可自訂，循環使用）
point_colors <- c("#E74C3C", "#9B59B6", "#F39C12", "#1ABC9C", "#3498DB")

# ============================================
# 讀取資料
# ============================================
df <- read_csv("20251212-0644.csv")

# Haversine 公式計算距離（公里）
haversine_distance <- function(lat1, lon1, lat2, lon2) {
  R <- 6371
  lat1_rad <- lat1 * pi / 180
  lat2_rad <- lat2 * pi / 180
  delta_lat <- (lat2 - lat1) * pi / 180
  delta_lon <- (lon2 - lon1) * pi / 180
  a <- sin(delta_lat / 2)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(delta_lon / 2)^2
  c <- 2 * atan2(sqrt(a), sqrt(1 - a))
  return(R * c)
}

# ============================================
# 計算每個測站到所有點位的距離
# ============================================
for (i in 1:nrow(target_points)) {
  col_name <- paste0("dist_", i)
  df[[col_name]] <- mapply(
    haversine_distance,
    target_points$lat[i], target_points$lon[i],
    df$latitude, df$longitude
  )
}

# 計算每個測站到最近點位的距離
dist_cols <- paste0("dist_", 1:nrow(target_points))
df$min_distance <- apply(df[, dist_cols], 1, min)
df$nearest_point <- apply(df[, dist_cols], 1, which.min)
df$nearest_point_name <- target_points$name[df$nearest_point]

# ============================================
# 分類測站（依據最近點位的距離）
# ============================================
df$category <- ifelse(
  df$min_distance <= radius_inner, 
  "inner",
  ifelse(
    df$min_distance <= radius_outer, 
    "middle",
    "outer"
  )
)

# ============================================
# 統計分析
# ============================================
stations_inner <- sum(df$category == "inner")
stations_middle <- sum(df$category == "middle")
stations_outer <- sum(df$category == "outer")

# 不重複站點統計（內圈 + 中圈）
unique_stations_in_range <- df %>% 
  filter(category %in% c("inner", "middle")) %>%
  distinct(station_id, .keep_all = TRUE)

unique_count_inner <- df %>% filter(category == "inner") %>% distinct(station_id) %>% nrow()
unique_count_middle <- df %>% filter(category == "middle") %>% distinct(station_id) %>% nrow()
unique_count_total <- nrow(unique_stations_in_range)

# 各點位統計
cat("\n")
cat("╔══════════════════════════════════════════════════════════════╗\n")
cat("║          空品測站距離篩選分析報告                            ║\n")
cat("╠══════════════════════════════════════════════════════════════╣\n")
cat("║ 參數設定                                                     ║\n")
cat("╠══════════════════════════════════════════════════════════════╣\n")
cat(sprintf("║  內圈半徑: %s 公里 (紅色)                                    ║\n", radius_inner))
cat(sprintf("║  外圈半徑: %s 公里 (藍色)                                    ║\n", radius_outer))
cat(sprintf("║  點位數量: %s 個                                             ║\n", nrow(target_points)))
cat("╠══════════════════════════════════════════════════════════════╣\n")
cat("║ 點位清單                                                     ║\n")
cat("╠══════════════════════════════════════════════════════════════╣\n")

for (i in 1:nrow(target_points)) {
  # 計算該點位範圍內的站點數
  inner_count <- sum(df[[paste0("dist_", i)]] <= radius_inner)
  outer_count <- sum(df[[paste0("dist_", i)]] <= radius_outer & df[[paste0("dist_", i)]] > radius_inner)
  cat(sprintf("║  %d. %s\n", i, target_points$name[i]))
  cat(sprintf("║     座標: (%.6f, %.6f)\n", target_points$lat[i], target_points$lon[i]))
  cat(sprintf("║     %dkm內: %d站 | %d-%dkm: %d站\n", 
              radius_inner, inner_count, radius_inner, radius_outer, outer_count))
  cat("║\n")
}

cat("╠══════════════════════════════════════════════════════════════╣\n")
cat("║ 📊 不重複站點統計（所有點位合併計算）                        ║\n")
cat("╠══════════════════════════════════════════════════════════════╣\n")
cat(sprintf("║  🔴 %dkm內不重複站點:     %4d 個                            ║\n", radius_inner, unique_count_inner))
cat(sprintf("║  🔵 %d-%dkm內不重複站點:   %4d 個                            ║\n", radius_inner, radius_outer, unique_count_middle))
cat(sprintf("║  ══════════════════════════                                 ║\n"))
cat(sprintf("║  ✅ 總計不重複站點:       %4d 個                            ║\n", unique_count_total))
cat(sprintf("║  🟢 範圍外站點:           %4d 個                            ║\n", stations_outer))
cat("╠══════════════════════════════════════════════════════════════╣\n")
cat(sprintf("║  📍 測站總數:             %4d 個                            ║\n", nrow(df)))
cat("╚══════════════════════════════════════════════════════════════╝\n")

# ============================================
# 建立 leaflet 地圖
# ============================================

# 計算地圖中心點（所有點位的平均）
center_lat <- mean(target_points$lat)
center_lon <- mean(target_points$lon)

map <- leaflet() %>%
  # 底圖
  addTiles(group = "OpenStreetMap") %>%
  addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "衛星影像") %>%
  
  # 設定視圖中心
  
  setView(lng = center_lon, lat = center_lat, zoom = 13)

# ============================================
# 為每個點位繪製圓圈和標記
# ============================================
for (i in 1:nrow(target_points)) {
  pt <- target_points[i, ]
  color <- point_colors[(i - 1) %% length(point_colors) + 1]
  
  # 計算該點位範圍內的站點數
  inner_count <- sum(df[[paste0("dist_", i)]] <= radius_inner)
  total_count <- sum(df[[paste0("dist_", i)]] <= radius_outer)
  
  # 繪製外圈
  map <- map %>%
    addCircles(
      lng = pt$lon,
      lat = pt$lat,
      radius = radius_outer * 1000,
      color = color,
      weight = 2,
      fillColor = color,
      fillOpacity = 0.03,
      dashArray = "5,5",
      popup = paste0(
        "<b>", pt$name, " - 外圈</b><br>",
        "半徑: ", radius_outer, " 公里<br>",
        "範圍內測站: ", total_count, " 個"
      ),
      group = "外圈範圍"
    )
  
  # 繪製內圈
  map <- map %>%
    addCircles(
      lng = pt$lon,
      lat = pt$lat,
      radius = radius_inner * 1000,
      color = color,
      weight = 3,
      fillColor = color,
      fillOpacity = 0.08,
      popup = paste0(
        "<b>", pt$name, " - 內圈</b><br>",
        "半徑: ", radius_inner, " 公里<br>",
        "範圍內測站: ", inner_count, " 個"
      ),
      group = "內圈範圍"
    )
  
  # 標記點位中心
  map <- map %>%
    addMarkers(
      lng = pt$lon,
      lat = pt$lat,
      popup = paste0(
        "<b>🎯 ", pt$name, "</b><br>",
        "緯度: ", round(pt$lat, 6), "<br>",
        "經度: ", round(pt$lon, 6), "<br>",
        "<hr style='margin:5px 0;'>",
        "🔴 ", radius_inner, "km內: ", inner_count, " 站<br>",
        "🔵 ", radius_inner, "-", radius_outer, "km: ", total_count - inner_count, " 站"
      ),
      icon = makeIcon(
        iconUrl = "https://raw.githubusercontent.com/pointhi/leaflet-color-markers/master/img/marker-icon-red.png",
        iconWidth = 25, iconHeight = 41,
        iconAnchorX = 12, iconAnchorY = 41
      ),
      group = "點位中心"
    )
}

# ============================================
# 繪製測站標記
# ============================================

# 內圈測站（紅色）
map <- map %>%
  addCircleMarkers(
    data = df[df$category == "inner", ],
    lng = ~longitude,
    lat = ~latitude,
    radius = 9,
    color = "#DC3545",
    fillColor = "#DC3545",
    fillOpacity = 0.9,
    weight = 2,
    popup = ~paste0(
      "<b>", thing_name, "</b><br>",
      "城市: ", city, "<br>",
      "鄉鎮: ", township, "<br>",
      "座標: (", round(latitude, 6), ", ", round(longitude, 6), ")<br>",
      "<hr style='margin:5px 0;'>",
      "<span style='color: #DC3545; font-weight: bold;'>",
      "🔴 最近點位: ", nearest_point_name, "<br>",
      "距離: ", round(min_distance, 3), " 公里</span>"
    ),
    group = paste0(radius_inner, "km內測站")
  )

# 中圈測站（藍色）
map <- map %>%
  addCircleMarkers(
    data = df[df$category == "middle", ],
    lng = ~longitude,
    lat = ~latitude,
    radius = 7,
    color = "#007BFF",
    fillColor = "#007BFF",
    fillOpacity = 0.8,
    weight = 2,
    popup = ~paste0(
      "<b>", thing_name, "</b><br>",
      "城市: ", city, "<br>",
      "鄉鎮: ", township, "<br>",
      "座標: (", round(latitude, 6), ", ", round(longitude, 6), ")<br>",
      "<hr style='margin:5px 0;'>",
      "<span style='color: #007BFF; font-weight: bold;'>",
      "🔵 最近點位: ", nearest_point_name, "<br>",
      "距離: ", round(min_distance, 3), " 公里</span>"
    ),
    group = paste0(radius_inner, "-", radius_outer, "km測站")
  )

# 外圈測站（綠色）
map <- map %>%
  addCircleMarkers(
    data = df[df$category == "outer", ],
    lng = ~longitude,
    lat = ~latitude,
    radius = 5,
    color = "#28A745",
    fillColor = "#28A745",
    fillOpacity = 0.6,
    weight = 1,
    popup = ~paste0(
      "<b>", thing_name, "</b><br>",
      "城市: ", city, "<br>",
      "鄉鎮: ", township, "<br>",
      "座標: (", round(latitude, 6), ", ", round(longitude, 6), ")<br>",
      "<hr style='margin:5px 0;'>",
      "<span style='color: #28A745;'>",
      "🟢 最近點位: ", nearest_point_name, "<br>",
      "距離: ", round(min_distance, 3), " 公里</span>"
    ),
    group = paste0(radius_outer, "km外測站")
  )

# ============================================
# 圖層控制與圖例
# ============================================
map <- map %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap", "CartoDB", "衛星影像"),
    overlayGroups = c(
      "內圈範圍", 
      "外圈範圍", 
      "點位中心", 
      paste0(radius_inner, "km內測站"),
      paste0(radius_inner, "-", radius_outer, "km測站"),
      paste0(radius_outer, "km外測站")
    ),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  
  # 圖例
  addLegend(
    position = "bottomright",
    colors = c("#DC3545", "#007BFF", "#28A745"),
    labels = c(
      paste0("🔴 ", radius_inner, "km內 (", unique_count_inner, " 不重複)"),
      paste0("🔵 ", radius_inner, "-", radius_outer, "km (", unique_count_middle, " 不重複)"),
      paste0("🟢 ", radius_outer, "km外 (", stations_outer, ")")
    ),
    title = paste0(
      "測站分布<br>",
      "<small>點位: ", nrow(target_points), " | ",
      "總站: ", nrow(df), "<br>",
      "<b>不重複: ", unique_count_total, "</b></small>"
    ),
    opacity = 0.9
  ) %>%
  
  # 比例尺
  addScaleBar(position = "bottomleft", options = scaleBarOptions(imperial = FALSE))

# ============================================
# 儲存地圖與結果
# ============================================
saveWidget(map, "stations_map_multipoint.html", selfcontained = TRUE)

# 輸出不重複站點清單
cat("\n📋 範圍內不重複站點清單:\n")
cat("─────────────────────────────────────────────────────────────\n")
unique_stations_in_range %>%
  select(station_id, thing_name, city, township, min_distance, nearest_point_name, category) %>%
  arrange(min_distance) %>%
  print(n = 50)

cat("\n✅ 地圖已儲存為 stations_map_multipoint.html\n")

🌟全文可以至下方連結觀看或是補充

全文分享至

https://www.facebook.com/LHB0222/

https://www.instagram.com/ahb0222/

有疑問想討論的都歡迎於下方留言

喜歡的幫我分享給所有的朋友 \o/

有所錯誤歡迎指教

# [:page_with_curl: 全部文章列表](https://hackmd.io/@LHB-0222/AllWritings)