基于sf：如何将公里级圆形叠加到SHP地图并匹配坐标系？

你手动用公式转换经纬度到公里的方式存在精度误差（地球是椭球体，不同纬度的缩放系数差异明显），更规范的做法是通过投影转换处理，同时可以直接用sf包生成匹配坐标系的圆形。以下分两种场景给出解决方案：

场景1：基于你已生成的plot.owin圆形，适配SHP坐标系 #

1. 先获取SHP图层的坐标系信息：

   library(sf)
   shp_layer <- st_read("你的文件路径.shp")
   shp_crs <- st_crs(shp_layer)
   

2. 将owin对象转换为sf格式，并补全原始坐标系信息（你的手动计算基于经纬度，先设为WGS84，EPSG:4326）：

   library(spatstat.geom)
   # 假设你已生成circle_owin对象
   circle_sp <- as(circle_owin, "SpatialPolygons")
   circle_sf <- st_as_sf(circle_sp)
   st_crs(circle_sf) <- 4326
   

3. 将圆形转换为SHP图层的坐标系：

   circle_sf_proj <- st_transform(circle_sf, shp_crs)
   

4. 叠加绘制：

   plot(st_geometry(shp_layer))
   plot(st_geometry(circle_sf_proj), add = TRUE, col = adjustcolor("blue", alpha.f = 0.3), lwd = 2)
   

场景2：用sf直接生成精准的10公里圆形（推荐） #

完全无需手动计算坐标，直接基于经纬度生成匹配坐标系的圆形：

1. 读取SHP并获取坐标系：

   shp_layer <- st_read("你的文件路径.shp")
   shp_crs <- st_crs(shp_layer)
   

2. 创建城市经纬度点并设置原始坐标系：

   # 替换为你的城市经纬度
   city_point <- st_sfc(st_point(c(城市经度, 城市纬度)), crs = 4326)
   

3. 转换到SHP的坐标系，生成10公里缓冲区（dist单位为米）：

   city_point_proj <- st_transform(city_point, shp_crs)
   circle_sf <- st_buffer(city_point_proj, dist = 10000)
   

4. 叠加绘制：

   plot(st_geometry(shp_layer), main = "城市10公里范围叠加")
   plot(st_geometry(circle_sf), add = TRUE, col = adjustcolor("red", alpha.f = 0.3), lwd = 2)
   

关键说明 #

• st_crs仅用于设置或读取坐标系信息，不能直接转换坐标。必须先给你的几何对象（点、圆形）指定原始坐标系，再用st_transform转换到目标坐标系。
• 手动坐标转换公式仅适用于粗略估算，高纬度地区误差会显著放大，建议始终使用投影转换工具保证精度。

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