时间戳

首页 > 解决方案 > 在 r/sf 中投影一个半径为 50 公里的四分之一圆？

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

问题描述

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

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

纽约州地图

我有两个问题：

  1. 我怎样才能有意义地设置这些圆的半径？有没有办法从投影的 CRS 中的坐标绘制一定距离（以公里为单位）的形状？到目前为止，我正在根据经度和纬度来定义半径，但距离会更有用。
  2. 在将它们投影并在 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")

标签： rmapsgeospatialsf

解决方案

对于任何对使用 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()

推荐阅读