在資料分析工作中，我們經常需要根據地理邊界來篩選資料。本文將完整示範如何下載鄉鎮邊界圖資，提取特定區域（以雲林縣斗六市為例），並篩選出落在該區域內的資料點。
前言：為什麼需要空間篩選？
在氣象資料、社會經濟調查、商業分析等領域，我們常常需要將資料點對應到特定的行政區域。例如：
- 找出特定城市內的氣象觀測站資料
- 分析某區域內的商業活動分布
- 統計鄉鎮層級的社會經濟指標
步驟一：環境準備
- 首先安裝及載入必要的 R 套件：
install.packages(c("sf", "dplyr", "ggplot2"))
library(sf) # 空間資料處理
library(dplyr) # 資料處理
library(ggplot2) # 資料視覺化
- 至政府資料開放平台，下載鄉鎮市區邊界圖資(網址：https://data.gov.tw/dataset/7441)
步驟二：讀取並篩選斗六市邊界
# 讀取鄉鎮市區界線檔案
town <- st_read("town_shapefile/TOWN_MOI_1120825.shp") %>%
# 篩選雲林縣斗六市
filter(COUNTYNAME == "雲林縣" & TOWNNAME == "斗六市") %>%
# 轉換座標系統為 WGS84 (經緯度)
st_transform(4326)
# 查看斗六市的基本資訊
print(paste("斗六市面積：", round(as.numeric(st_area(town))/1000000, 2), "平方公里"))
print(paste("村里數：", nrow(town)))
步驟三：準備範例資料集
假設我們有一個包含經緯度和溫度的資料集：
# 建立範例資料（實際應用中這可能是你的觀測資料）
set.seed(123)
data <- tibble(
year = rep(2023, 1000),
LON = runif(1000, 120.4, 120.7), # 斗六市周邊經度範圍
LAT = runif(1000, 23.6, 23.8), # 斗六市周邊緯度範圍
年平均溫度 = runif(1000, 22, 28)
)
print(paste("原始資料筆數：", nrow(data)))
步驟四：將資料轉換為空間物件
# 將普通資料框轉換為空間資料
data_sf <- st_as_sf(data,
coords = c("LON", "LAT"), # 指定經緯度欄位
crs = 4326) # 設定座標系統
# 確認座標系統一致
st_crs(data_sf) == st_crs(town)
步驟五：進行空間篩選
有幾種方法可以篩選落在斗六市內的點：
方法1：使用 st_join
# st_join() 會找出所有落在 town 多邊形範圍內的點位
# left = FALSE 確保只保留落在範圍內的點
data1 <- st_join(data_sf, town, left = FALSE)
方法2：使用 st_intersects
# 找出與斗六市邊界相交的點
data2 <- st_intersects(data_sf, town, sparse = FALSE)
data2 <- data_sf[data2, ]
print(paste("斗六市內的資料筆數：", nrow(data2)))
步驟六：視覺化驗證結果
# 建立地圖驗證篩選結果
ggplot() +
# 繪製斗六市邊界
geom_sf(data = town, fill = "lightblue", alpha = 0.5) +
# 繪製所有原始資料點
geom_sf(data = data_sf, color = "gray", size = 1, alpha = 0.5) +
# 繪製斗六市內的點
geom_sf(data = data1, color = "red", size = 1) +
# 設定標題
labs(title = "斗六市範圍內資料點篩選結果",
subtitle = paste("總點數：", nrow(data_sf), "，斗六市內：", nrow(douliu_data))) +
theme_minimal()
步驟七：轉回一般資料框進行後續分析
# 將空間資料轉回普通資料框
data1_df <- data1 %>%
st_drop_geometry() # 移除空間資訊
# 現在可以進行一般的資料分析了
summary(data1_df$年平均溫度)
# 計算斗六市內的年平均溫度統計
data1_summary <- data1_df %>%
group_by(year) %>%
summarise(
平均溫度 = mean(年平均溫度),
最高溫度 = max(年平均溫度),
最低溫度 = min(年平均溫度),
資料點數 = n()
)
print(data1_summary)
實際應用案例
案例1：氣象資料分析
# 篩選特定年份的氣象觀測站資料
weather_stations <- read.csv("氣象站資料.csv") %>%
st_as_sf(coords = c("經度", "緯度"), crs = 4326)
douliu_stations <- weather_stations[st_intersects(weather_stations, town, sparse = FALSE), ]
案例2：商業據點分析
# 分析斗六市內的便利商店分布
convenience_stores <- read.csv("便利商店位置.csv") %>%
st_as_sf(coords = c("LON", "LAT"), crs = 4326)
douliu_stores <- convenience_stores[st_intersects(convenience_stores, town, sparse = FALSE), ]
常見問題與解決方案
Q1：座標系統不一致怎麼辦？
# 確認兩者的座標系統
st_crs(data_sf)
st_crs(town)
# 如果不一致，進行轉換
town <- st_transform(town, st_crs(data_sf))
Q2：如何處理跨邊界的點？
# 如果需要精確計算點在區域內的面積比例
area_ratio <- st_intersection(data_sf, town) %>%
mutate(area_ratio = as.numeric(st_area(geometry)) / as.numeric(st_area(data_sf)))
