R语言根据经纬度获得对应样本的省份

一.输入

包含样本的经纬度

二.代码

library(sf)
library(dplyr)
library(osmextract)# 1. 准备省份中英文对照表
province_names <- data.frame(name_en = c("Beijing", "Tianjin", "Shanghai", "Chongqing", "Hebei", "Shanxi", "Inner Mongolia", "Liaoning","Jilin", "Heilongjiang", "Jiangsu", "Zhejiang","Anhui", "Fujian", "Jiangxi", "Shandong","Henan", "Hubei", "Hunan", "Guangdong","Guangxi", "Hainan", "Sichuan", "Guizhou","Yunnan", "Tibet", "Shaanxi", "Gansu","Qinghai", "Ningxia", "Xinjiang", "Taiwan"),province_cn = c("北京市", "天津市", "上海市", "重庆市","河北省", "山西省", "内蒙古自治区", "辽宁省","吉林省", "黑龙江省", "江苏省", "浙江省","安徽省", "福建省", "江西省", "山东省","河南省", "湖北省", "湖南省", "广东省","广西壮族自治区", "海南省", "四川省", "贵州省","云南省", "西藏自治区", "陕西省", "甘肃省","青海省", "宁夏回族自治区", "新疆维吾尔自治区", "台湾省")
)# 2. 读取元数据
metadata <- read.csv("lib_metadata2.csv") %>% filter(!is.na(latitude), !is.na(longitude))# 3. 获取中国行政区划数据（两种方式任选其一）
# 方式A：使用osmextract获取数据
poly_china <- openstreetmap_fr_zones %>% filter(parent == "china") %>% left_join(province_names, by = c("name" = "name_en")) %>% mutate(province_cn = ifelse(is.na(province_cn), name, province_cn))# 方式B：使用本地GeoJSON文件（推荐，更稳定）
# china_map <- st_read("D:/path/to/china.json") %>% 
#   left_join(province_names, by = c("name" = "name_en")) %>% 
#   mutate(province_cn = coalesce(province_cn, name))# 4. 转换为空间点数据
points <- st_as_sf(metadata,coords = c("longitude", "latitude"),crs = 4326,  # WGS84坐标系remove = FALSE  # 保留原始列
) %>% st_transform(st_crs(poly_china))  # 转换为与多边形相同的坐标系# 5. 执行空间连接（带中文省份名称）
result <- st_join(points, poly_china, join = st_within) %>% select(lib_id, genus, species, latitude, longitude, region_en = name, region_cn = province_cn)# 6. 处理未匹配的点（可选）
if(any(is.na(result$region_cn))) {# 方法A：标记未匹配点result <- result %>% mutate(region_cn = ifelse(is.na(region_cn), "未知地区", region_cn))# 方法B：使用最近邻匹配（更精确）# unmatched <- result %>% filter(is.na(region_cn))# nearest <- st_nearest_feature(unmatched, poly_china)# result[is.na(result$region_cn), "region_cn"] <- poly_china$province_cn[nearest]
}# 7. 检查结果
print(result %>% count(region_cn, sort = TRUE) %>% as.data.frame())# 8. 保存结果（UTF-8编码支持中文）
write.csv(result, "D:/file/BGI/蚊虫项目/Wolbachia/bam/metadata_with_regions.csv", row.names = FALSE, fileEncoding = "UTF-8")

结果