如何用R将县坐标匹配到Shapefile多边形并实现地图着色？

嘿，我来帮你搞定这个问题！你现在需要把县的经纬度数据匹配到Shapefile的多边形里，然后给不同县对应的区域上色对吧？下面我给你两种方案，一种基于你现在用的工具，另一种是更现代的推荐做法，方便你选择：

1. 基于你现有代码的解决方案（maptools + sp） #

你目前用的maptools虽然可以实现，但需要补充空间匹配的步骤，具体操作如下：

步骤1：将县数据转为空间点对象并匹配多边形 #

首先要把你的县列表转换成空间点格式，再用over()函数判断每个县的点属于哪个多边形：

library(maptools)
library(sp)
library(ggplot2)

# 读取你的Shapefile和县数据
area <- readShapePoly("Path/file.shp")
counties <- read.csv("./Path/counties.csv")

# 将县数据转为空间点对象，关键是要和Shapefile用**相同的坐标系**！
coordinates(counties) <- ~long + lat
proj4string(counties) <- proj4string(area)  # 同步坐标系，避免匹配错误

# 用over()函数匹配每个县所属的多边形，这里的`ID`是Shapefile多边形的标识列，你要换成实际列名（比如NAME、COUNTY_ID等）
counties$poly_id <- over(counties, area)$ID

步骤2：合并数据并绘制带颜色的地图 #

接下来把多边形的坐标数据和匹配后的县信息关联，就能给不同县的区域上色了：

# 把Shapefile的属性表转成数据框
area_df <- as.data.frame(area)
# 用匹配的poly_id关联多边形属性和县数据
area_df <- merge(area_df, counties, by.x = "ID", by.y = "poly_id")

# 把多边形转为ggplot可用的格式，同时保留多边形ID
area.points <- fortify(area, region = "ID")
# 合并坐标数据和属性数据
map_data <- merge(area.points, area_df, by.x = "id", by.y = "ID")

# 绘制地图，这里用`county_name`作为上色依据，你可以换成自己的县标识列
map <- ggplot() +
  geom_polygon(data = map_data, aes(x = long, y = lat, group = group, fill = county_name), color = 'gray', size = .2) +
  coord_map()  # 优化地图投影，让形状更准确
print(map)

2. 推荐的现代方案（sf包） #

maptools已经停止维护了，现在处理空间数据的标准工具是sf包，代码更简洁且不易出错：

library(sf)
library(ggplot2)

# 读取Shapefile和县数据，sf包直接读取就是空间格式
area_sf <- st_read("Path/file.shp")
counties <- read.csv("./Path/counties.csv")

# 将县数据转为空间点对象，自动同步坐标系
counties_sf <- st_as_sf(counties, coords = c("long", "lat"), crs = st_crs(area_sf))

# 匹配每个县所属的多边形，st_join会直接把多边形的属性关联到点上
counties_with_poly <- st_join(counties_sf, area_sf)

# 把县的信息合并到多边形数据里，假设用`county_name`作为关联字段
area_with_counties <- merge(area_sf, counties_with_poly, by = "county_name")

# 绘制地图，sf和ggplot完美兼容，直接用geom_sf即可
map <- ggplot() +
  geom_sf(data = area_with_counties, aes(fill = county_name), color = 'gray', size = .2) +
  theme_minimal()
print(map)

重要注意事项 #

• 务必保证坐标系一致：如果Shapefile和县数据用的是不同的坐标系（比如一个是WGS84，一个是UTM），匹配结果会完全错误，用proj4string()（sp包）或st_crs()（sf包）查看并同步坐标系。
• 确认多边形的标识列：你需要知道Shapefile里用来区分不同多边形的字段（比如NAME、ID），不要用错列名导致匹配失败。
• 处理边界点情况：如果某个县的点正好落在多边形边界上，可能会匹配到多个多边形，这时可以用st_intersection（sf包）做更精确的空间判断，或者取第一个匹配结果。

内容的提问来源于stack exchange，提问作者Peter