基于经纬度识别一级行政单元并可视化标红的技术实现求助

需求 #

从数据集的经纬度信息中识别对应的一级行政单元，并在可视化地图中将该行政单元填充为红色。

原错误代码 #

library(rnaturalearth)
library(sf)
library(dplyr)
library(ggplot2)


find_admin_areas <- function(lon, lat) {
  # 
  all_states <- ne_states(returnclass = "sf")
  
  # 
  point_df <- st_as_sf(data.frame(lon = lon, lat = lat), coords = c("lon", "lat"), crs = st_crs(all_states))
  
  # 
  admin_areas <- st_join(point_df, all_states)
  
  # 
  admin_areas <- st_make_valid(admin_areas)
  
  return(admin_areas)
}

# 
locations_df <- data.frame(
  location = c("Kampong Chhnang", "Prey Veng", "Bokeo", "Xieng Khouang", "Phatthalung", "Chachoengsao",
               "Kamphaeng Phet", "Maha Sarakham", "Nonthaburi", "Phitsanulok", "Roi Et", "Si Sa Ket",
               "Surin", "Uthai Thani", "Lao Cai", "Thua Thien-Hue", "Phang Nga", "An Giang"),
  lon = c(104.5598351, 105.4249716, 101.8603, 103.1700, 100.0694874, 101.4314805,
          99.53470511, 103.1683362, 100.394886, 100.5448427, 103.8151837, 104.3711179,
          103.658002, 99.47934782, 103.9768, 107.501161, 98.42208133, 105.182631),
  lat = c(12.16634032, 11.39794884, 20.2672, 19.4326, 7.510663288, 13.60579712,
          16.33080535, 15.9978543, 13.92069527, 16.98237825, 15.91654777, 14.85512615,
          14.88487277, 15.34854839, 22.3820, 16.32788631, 8.4500, 10.5143)
)

# 
admin_areas <- locations_df %>%
  rowwise() %>%
  mutate(admin_area = list(find_admin_areas(lon, lat))) %>%
  select(location, admin_area)
## <- Error code 

# 
ggplot() +
  geom_sf(data = ne_states(), fill = "lightgray", color = "white") +  
  geom_sf(data = admin_areas$admin_area, fill = "red", color = "red", alpha = 0.5) + 
  labs(title = "Administrative Areas in Southeast Asia") +  
  theme_minimal()  

错误原因与修复方案 #

1. 重复加载数据集：原函数每次调用都重新加载ne_states()，既浪费资源又容易出错，应提前加载并复用。
2. 空间连接返回结果错误：st_join返回的是带行政属性的点数据，而非行政单元的面几何，需要根据匹配到的唯一标识（如adm1_code）从ne_states中提取对应面。
3. 数据结构不兼容ggplot：rowwise生成的嵌套数据框无法被geom_sf直接识别，需扁平化处理为标准sf对象。
4. CRS显式指定：创建点时显式指定EPSG:4326，避免坐标系匹配歧义。

修正后的完整代码 #

library(rnaturalearth)
library(sf)
library(dplyr)
library(ggplot2)
library(tidyr)

# 提前加载全球一级行政单元数据，只加载一次
all_states <- ne_states(returnclass = "sf") %>%
  st_make_valid()  # 提前处理几何有效性

# 定义匹配函数：根据经纬度返回对应的行政单元面
find_admin_area <- function(lon, lat) {
  # 创建点并指定坐标系为WGS84
  point <- st_as_sf(data.frame(lon = lon, lat = lat), 
                    coords = c("lon", "lat"), 
                    crs = 4326)
  
  # 空间连接匹配行政单元属性
  matched_attr <- st_join(point, all_states, join = st_within)
  
  # 根据匹配到的adm1_code提取对应的行政单元面
  if (!is.na(matched_attr$adm1_code)) {
    admin_area <- all_states %>%
      filter(adm1_code == matched_attr$adm1_code)
  } else {
    admin_area <- NULL
  }
  
  return(admin_area)
}

# 原始位置数据
locations_df <- data.frame(
  location = c("Kampong Chhnang", "Prey Veng", "Bokeo", "Xieng Khouang", "Phatthalung", "Chachoengsao",
               "Kamphaeng Phet", "Maha Sarakham", "Nonthaburi", "Phitsanulok", "Roi Et", "Si Sa Ket",
               "Surin", "Uthai Thani", "Lao Cai", "Thua Thien-Hue", "Phang Nga", "An Giang"),
  lon = c(104.5598351, 105.4249716, 101.8603, 103.1700, 100.0694874, 101.4314805,
          99.53470511, 103.1683362, 100.394886, 100.5448427, 103.8151837, 104.3711179,
          103.658002, 99.47934782, 103.9768, 107.501161, 98.42208133, 105.182631),
  lat = c(12.16634032, 11.39794884, 20.2672, 19.4326, 7.510663288, 13.60579712,
          16.33080535, 15.9978543, 13.92069527, 16.98237825, 15.91654777, 14.85512615,
          14.88487277, 15.34854839, 22.3820, 16.32788631, 8.4500, 10.5143)
)

# 批量匹配并合并结果
matched_admin_areas <- locations_df %>%
  rowwise() %>%
  mutate(admin_area = list(find_admin_area(lon, lat))) %>%
  unnest(admin_area) %>%  # 扁平化嵌套数据
  st_as_sf()  # 转换为sf对象

# 可视化
ggplot() +
  geom_sf(data = all_states, fill = "lightgray", color = "white", linewidth = 0.2) +  
  geom_sf(data = matched_admin_areas, fill = "red", color = "darkred", alpha = 0.6) + 
  labs(title = "东南亚一级行政单元匹配结果") +  
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

内容的提问来源于stack exchange，提问作者新子泰平