在 r/sf 中投影一个半径为 50 公里的四分之一圆？答案

【问题标题】：Projecting a quartered circle with a 50km radius in r/sf?在 r/sf 中投影一个半径为 50 公里的四分之一圆？
【发布时间】：2021-09-30 01:46:55
【问题描述】：

我希望创建一系列四等分的圆（即分成 4 个相等的象限），每个圆的半径为 50 公里，我可以将它们映射到整个美国的各种经度和纬度。我还想要根据需要旋转这些四等圆的选项。

使用下面的代码（以及来自here 的指导），我已经能够进行以下操作：

New York State Map

我有两个问题：

如何有意义地设置这些圆的半径？有没有办法从投影的 CRS 中的坐标绘制一定距离（以公里为单位）的形状？到目前为止，我根据经度和纬度来定义半径，但距离会更有用。
在 WGS84 中投影和映射后，我的圆圈似乎变成了椭圆。有什么办法可以防止这种情况发生吗？

我很乐意考虑其他方法。谢谢！

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

#Two functions to create coordinate quartered circle polygons

#x = long, y = lay, r = radius, theta_rotate = rotation
st_wedge <- function(x,y,r,start,width, theta_rotate){
   n <- 20
      theta = seq(start+theta_rotate, start+width+theta_rotate, 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(x, y, r, nsegs, theta_rotatex){
   width = (2*pi)/nsegs
   starts = (1:nsegs)*width
   polys = lapply(starts, function(s){st_wedge(x,y,r,s,width, theta_rotatex)})
   #Cast to crs 4326, WGS84
   mpoly = st_cast((st_sfc(polys, crs = 4326)), "MULTIPOLYGON")
   mpoly
}

#Create quartered sf circle polygon 
custom_circle_sf <- st_wedges(x = -76, y = 43, r = .3, nsegs = 4, theta_rotatex = 200) %>%
 st_sf() %>% 
   mutate(group = row_number()) %>% dplyr::select(group, geometry)

#Create New York State sf polygon
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")

#Plot results
ggplot() +
   geom_sf(data=ny_map_sf,
           size = 1,
           colour = "blue", 
           fill   = "white") + 
   geom_sf(data=custom_circle_sf,
           size = .1,
           aes(fill=group),
           colour = "white")

【问题讨论】：

我还没有运行你的代码，但看起来你让它变得比它需要的更复杂。您可以将您的几何图形转换为以米为单位的美国 CRS —EPSG:4269 是一。然后只需使用sf::st_buffer 制作一个大小合适的圆圈。可能是一种更简单的方法，也可以使用一些sf 函数将其分解为楔形，可能是通过在缓冲区和边界框部分之间进行交叉。我理解对了吗？
另外，最好将第二个问题分成自己的帖子
谢谢卡米尔。我被认为是沿着这些思路的解决方案（您理解正确）。我相信您所描述的方法的想法是 1）创建一个以圆心为中心的 2x2“棋盘”； 2）给每个复选框一个等于圆半径的长度；和 3) 找到一种方法，使这个棋盘围绕原点（即圆心）旋转给定的度数。
我想我可以使用 st_bbox 创建边界框，然后使用 st_split 将框分成四个象限。还有什么听起来更可取的吗？
是的，除了使用st_buffer 来获得圈子，这就是我的想象。不确定我能否找到完整的答案，但我会继续考虑这个问题。

标签： r maps geospatial sf

【解决方案1】：

对于任何对使用 R 在 sf 中分割多边形感到好奇的人，这就是我解决这个问题的方法：

#Function to create circle with quadrants. Save desired projection as projected_crs
create_circle <- function(lat_x, long_y, theta_x, buffer_m){
   #Create circle with radius buffer_m centered at (lat_x, long_y)
   circle_buffer  <- st_point(c(lat_x, long_y)) %>%  st_sfc(crs = 4326) %>% 
      st_cast("POINT")  %>% 
      st_transform(projected_crs) %>%
      st_buffer(buffer_m)
   
   #Create two orthogonal lines at origin 
   p1 <- rbind(c(lat_x,long_y - 1), c(lat_x,long_y + 1))
   p2 <- rbind(c(lat_x+1,long_y), c(lat_x-1,long_y))
   mls <- st_multilinestring(list(p1,p2))  %>% st_sfc(crs = 4326) %>% 
      st_transform(projected_crs) 
   
   #Use orthogonal lines to split circle into 4 quadrants
   x1 <- st_split(circle_buffer, mls) 
   
   #Convert origin into projected CRS
   center_in_crs  <- st_point(c(lat_x, long_y)) %>% 
      st_sfc(crs = 4326) %>%
      st_transform(projected_crs)
   
   sp_obj <- x1 %>% st_collection_extract(type="POLYGON") %>%
      #Convert to spatial to use sp functions
      as_Spatial() %>% 
      #rotate x degrees
      elide(rotate = theta_x + 45, center  = center_in_crs[[1]]) %>% 
      #return to sf 
      st_as_sf()

【讨论】：

【解决方案2】：

关于您的问题 2：“圆圈似乎变成了椭圆”。如果将 coord_equal() 函数添加到 ggplot 中，则网格将是正方形，而椭圆将显示为圆形。

【讨论】：

谢谢！这很有帮助