R sf包绘制50km半径可旋转四等分圆及投影变形问题咨询

核心思路 #

你的两个问题本质都是在地理坐标系（WGS84/EPSG:4326）下做平面几何计算导致的，统一切换到以米为单位的北美投影坐标系（EPSG:5070）即可同时解决：

• 该投影覆盖美国全境，距离计算误差小，默认单位为米，可直接传入50000对应50km半径
• 投影后平面几何是正形的，生成的圆不会出现变形

问题对应说明 #

1. 固定公里半径设置方法
   不要用经纬度度数定义半径，1度经度对应的实际距离随纬度升高不断缩小，北纬40度附近1度经度仅对应约85km，误差极大。正确流程是：

• 先将中心经纬度点转换为投影坐标系
• 直接用米为单位设置半径做几何计算
• 最终需要WGS84坐标的话，计算完成后再转换回EPSG:4326即可

2. WGS84下椭圆变形解决方法
   WGS84是球面坐标，平面展示时默认经纬度等比例映射，高纬度区域东西方向被压缩，天然会把圆拉成椭圆。只要生成几何、绘制环节都使用同一投影坐标系，即可完全避免变形。

修正后完整代码 #

library(sf)
library(ggplot2)
library(maps)
library(dplyr)

# 修正后的扇形生成函数：x,y为投影坐标系下的坐标，r单位为米，theta_rotate单位为角度
st_wedge <- function(x,y,r,start,width, theta_rotate){
  n <- 20
  # 角度转弧度适配三角函数计算
  theta_rotate_rad <- theta_rotate * pi / 180
  theta = seq(start+theta_rotate_rad, start+width+theta_rotate_rad, length=n)
  xarc = x + r*sin(theta) 
  yarc = y + r*cos(theta)
  xc = c(x, xarc, x) 
  yc = c(y, yarc, y)
  st_polygon(list(cbind(xc,yc)))   
}

st_wedges <- function(lon, lat, r_m, nsegs, theta_rotatex){
  # 先创建WGS84下的中心点，再转换为北美等面积投影EPSG:5070（单位：米）
  center_pt <- st_sfc(st_point(c(lon, lat)), crs = 4326) %>% 
    st_transform(5070)
  center_coord <- st_coordinates(center_pt)
  x <- center_coord[1]
  y <- center_coord[2]
  
  width = (2*pi)/nsegs
  starts = (1:nsegs)*width
  polys = lapply(starts, function(s){st_wedge(x,y,r_m,s,width, theta_rotatex)})
  # 生成投影坐标系下的多面，可按需转回到WGS84
  mpoly = st_sfc(polys, crs = 5070) %>% 
    st_cast("MULTIPOLYGON")
  # 如果需要输出WGS84坐标，取消下一行注释即可
  # mpoly <- st_transform(mpoly, 4326)
  mpoly
}

# 生成50km半径的四等分圆，传入中心经纬度、半径50000米、旋转角度200
custom_circle_sf <- st_wedges(lon = -76, lat = 43, r_m = 50000, nsegs = 4, theta_rotatex = 200) %>%
  st_sf() %>% 
  mutate(group = row_number()) %>% 
  dplyr::select(group, geometry)

# 纽约州地图也转换为同一投影，避免绘制变形
ny_map_sf <- map_data("state", region="new york")  %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
  group_by(group) %>%
  summarise(geometry = st_combine(geometry)) %>% 
  st_cast("POLYGON") %>% 
  st_transform(5070)

# 绘制结果，完全是正圆形，半径也符合50km要求
ggplot() +
  geom_sf(data=ny_map_sf,
          size = 1,
          colour = "blue", 
          fill   = "white") + 
  geom_sf(data=custom_circle_sf,
          size = .1,
          aes(fill=as.factor(group)),
          colour = "white") +
  labs(fill = "象限")

其他优化建议 #

如果只需要做单个小区域的圆，也可以使用对应区域的UTM投影，距离精度会比EPSG:5070更高。如果需要批量生成美国全境的多个四等分圆，直接复用以上代码传入不同经纬度即可。

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