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

一种方法，使用 {sf} 和 {dplyr}：

编辑 根据 @Chris 和 @Salvador 的有用评论，此解决方案已大幅重写。

TLDR：下面的解决方案创建一个水平基线，缓冲此基线以产生所选宽度的水平条带，创建该条带的相邻克隆，将它们组合成一个多多边形，然后旋转和移动多个多边形以从所选锚向南扩展。

1. helper 函数返回相邻条带的多多边形：get_strips(...)

get_strips <- function(n = 10, ## number of adjacent strips
                       strip_width = 1e3, # unit: m
                       strip_length = 1e5,
                       crs = 3857 ## default: 'google' mercator
                       ){
  ## make horizontal line of length 'strip_length':
  baseline <- matrix(c(-strip_length/2, strip_length/2, 0, 0), ncol = 2) |>
    st_linestring()
  ## make line a polygon of width 'strip_width'
  basestrip  <- baseline |> 
    st_buffer(strip_width/2, endCapStyle = 'FLAT')
  ## return a spatial dataframe of n adjacent strips from north to south:
  st_sf(strip_id = sprintf('strip_%02.f', 1:n),
        geometry = Map(1:n,
                       f = \(index) basestrip - c(0, index * strip_width - strip_width/2)
                       )|>
          st_sfc() |> st_cast('MULTIPOLYGON') |> st_set_crs(crs)
        )
}

2. 用于旋转空间数据帧几何图形的帮助程序函数

    rotate_feature <- function(feature,
                               angle = 0, ## rotation in degree
                               anchor = c(0, 0) ## top center point of grid
                               ){
      crs <- st_crs(feature)
      a <- -pi * angle/180 ## convert to radiant ccw
      ## create rotation matrix:
      rotmat <- matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
      ## rotate, recenter, restore CRS:
      st_geometry(feature) <- st_geometry(feature) * rotmat + anchor
      feature |> st_set_crs(crs)
    }

3. 加载库，创建播放数据

    library(sf)
    library(dplyr)
    
    the_lake <- readRDS(file.choose()) ## spatial feature supplied by OP
    the_anchor <- c(5.08e5, 4.224e6) ## anchor strip grid here

4. 生成所需的空间要素

    the_strips <- 
      get_strips(crs = st_crs(the_lake)) |>
      rotate_feature(angle = 30, anchor = the_anchor)
    
    the_cropped_strips <- the_strips |>
      st_intersection(the_lake)
    
    ## casting polygons to multilinestrings to only
    ## obtain the corner points (not areas) on intersection:
    the_corners <- 
      st_intersection(
        st_cast(the_strips, 'MULTILINESTRING'),
        st_cast(the_lake, 'MULTILINESTRING')
      )

5. 检查

    library(ggplot2)
    
    anchor_point <- the_anchor |> st_point() |>
                                  st_sfc(crs = st_crs(the_lake))
    
    ggplot() +
      geom_sf(data = the_lake) +
      geom_sf(data = the_strips, fill = NA) +
      geom_sf(data = the_cropped_strips, aes(fill = strip_id)) +
      geom_sf(data = the_corners) +
      geom_sf(data = anchor_point, pch = '+', size = 5) +
      geom_sf_label(data = anchor_point, nudge_y = 1e3, label = 'Anchor') +
      coord_sf(xlim = 1e5 * c(5, 5.2), ylim = 1e6 * c(4.21, 4.235))

6. 提取每个多边形的角：

    the_corners |>
      rowwise() |>
      reframe(strip_id,
              coords = st_coordinates(geometry)
              )

    ## # A tibble: 50 x 2
    ##    strip_id coords[,1]     [,2]  [,3]
    ##    <chr>         <dbl>    <dbl> <dbl>
    ##  1 strip_01    508864. 4223922.     1
    ##  2 strip_01    507927. 4223381.     1
    ##  3 strip_01    508230. 4222400.     1
    ##  4 strip_01    509328. 4223035.     1
    ##  5 strip_02    508230. 4222400.     1
    ##  6 strip_02    508791. 4221570.     1

另一种使用空间填充多边形的方法，这里是正方形。

n100 = c(100,100)
bay_square_100 = st_make_grid(
bay,
cellsize = c(diff(st_bbox(bay)[c(1,3)]), diff(st_bbox(bay)[c(2,4)]))/n100,
offset = st_bbox(bay)[c('xmin', 'ymin')],
what = 'polygon',
square = TRUE
)
# get index of just the ones inside
bay_sq_100_cont = st_contains_properly(bay, bay_square_100)
contained_bay_sq = bay_square_100[unlist(bay_sq_100_cont)]

索引与非常有用的 ， 序列的 rle， 1,2,3...icontained_bay_sqcgwtools::seqle

end = cumsum(cgwtools::seqle(unlist(bay_sq_100_cont))$lengths)
start = end - cgwtools::seqle(unlist(bay_sq_100_cont))$lengths +1

这种方法会创建一个对象，直到你遇到 knapply 注释、提醒或通知有参数的对象之前，这个对象似乎是行不通的。st_unionunlistrecursive

unioned_rows = list()
for (i in 1:length(start)) {
unioned_rows[[i]] = st_union(c(bay_cont[start[i]:end[i]]))
}
unioned_row4 = st_sfc(unlist(unioned_rows, recursive = FALSE)) #many difficult objects created

plot(unioned_row4, col = c('red', 'blue', 'green', 'orange', 'pink') )

这可以通过六边形和沿对角线的索引来完成，以创建人造角度，但我认为@I_O的方法更优雅灵活。

实际上，网格可以旋转，如所示 Stackoverflow - 旋转网格

library(sf)
Linking to GEOS 3.12.0, GDAL 3.8.0, PROJ 9.3.0; sf_use_s2() is TRUE
 bay = st_geometry(readRDS('~/Downloads/bay.RDS'))
 rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
 tran = function(geo, ang, center) (geo - center) * rot(ang * pi / 180) + center
 center = st_centroid(st_union(bay))
 rotang = -30 # sign, -CounterClockWise +CW
 grd = st_make_grid(tran(bay, -rotang, center), cellsize = 100, crs = st_crs(bay))
 grd_rot = tran(grd, rotang, center)
 bay_rot_100_cont_idx = st_contains_properly(bay, grd_rot)
Error in st_geos_binop("contains_properly", x, y, sparse = sparse, prepared = TRUE,  : 
  st_crs(x) == st_crs(y) is not TRUE
 st_crs(grd_rot) = 26910
 bay_rot_100_cont_idx = st_contains_properly(bay, grd_rot)
 cont_bay_rot = grd_rot[unlist(bay_rot_100_cont_idx)]
 end_rot = cumsum(cgwtools::seqle(unlist(bay_rot_100_cont_idx))$lengths)
 start_rot = end_rot - cgwtools::seqle(unlist(bay_rot_100_cont_idx))$lengths + 1
 unioned_rows_rot = list()
 for (i in 1:length(start_rot)) {
 unioned_rows_rot[[i]] = st_union(c(cont_bay_rot[start_rot[i]:end_rot[i]]))
 }
 unioned_rot_sfc = st_sfc(unlist(unioned_rows_rot, recursive = FALSE)
plot(unioned_rot_sfc, col = c('green', 'white', 'red'))