如何从shapefile制作小多边形并提取坐标

Question

我这里有一个小形状文件：https://login.filesanywhere.com/fs/v.aspx?v=8c6d62865c626fb4a2ab称为bay.RDS

 library(tmap)
 library(leaflet)
 library(mapview)
 bay <- readRDS('bay.RDS')
 mapview(bay)

我正在尝试在下面的点从海岸到海岸绘制 50 米宽的多边形，并从中提取坐标。 我需要多个彼此相邻但独立绘制的多边形，因为我需要根据一列值对它们进行颜色编码

pt <- st_sfc(st_point(c(-122.91, 38.17)), crs=4326)
mapview(bay) + mapview(pt, color ='red')

我在下面添加了我想要完成的任务的快照。如果有人能告诉我如何做，我将不胜感激 制作几个多边形并从中提取坐标，以便我可以将它们合并到我的数据集中。如果不清楚请告诉我

Answer 1

一种方法，使用 {sf} 和一些三角函数：

• 创建绘制条带的点：

library(sf)

the_point <- 
  st_point(c(5.08e5, 4.224e6)) |>
  st_sfc(crs = st_crs(the_lake))

• 创建条带的辅助函数（仅限平行线，不是多边形）：

get_parallels <- function(point,
                          n = 5, ## number of strips (strip borders actually)
                          angle = 0, ## strip angle in degrees
                          strip_width = 1e3, strip_length = 10e3, 
                          crs ## CRS of the feature to intersect with
                          ){
  n = ceiling(n/2)
  a <- -pi * angle/180
  rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
  centroid <- point |> st_coordinates()
  endpoints <- c(st_point(c(-strip_length/2, 0)), st_point(c(strip_length/2, 0))) |>
    st_multipoint()
  baseline <- st_linestring(endpoints * rot(a) + centroid) ## rotated baseline through sample point
  ## create n parallels to baseline with normal distance strip_width
  Map(-n:n * strip_width, f = \(x) baseline + x * c(sin(a), sqrt(sin(a)^2 + cos(a)^2))) |>
    st_multilinestring() |>
    st_sfc(crs = crs)
}

• 创建条带，

  the_lake

  是通过

  sf

  您的数据获得的

  readRDS

  ：

the_parallels <- get_parallels(the_point, angle = 30, strip_width = 1000, crs = st_crs(the_lake))

• 测试图（未显示）

plot(the_lake |> st_geometry())
plot(the_point, add = TRUE)
plot(the_parallels, add = TRUE)

• 添加纬线和海岸线的交点（点）：

st_intersection(the_parallels,
                st_cast(the_lake, 'MULTILINESTRING')
                ) |>
  plot(add = TRUE, col = 'red')

结果：

编辑

要获取每个海岸到海岸轨迹的端点，您可以继续：

st_intersection(the_lake, the_parallels) |>
  st_cast('MULTIPOINT') |> 
  st_coordinates()

返回（输出被截断）：

             X       Y L1
 [1,] 509124.8 4220783  1
 [2,] 510599.9 4221635  1
 [3,] 508570.0 4221752  2
 [4,] 509700.3 4222404  2
## ...

...其中 L1 表示轨迹（两行 = 每个交叉点的端点）