如何设置锚点并从 shapefile 创建多边形

Question

我正在跟进这里的示例：如何从 shapefile 制作小多边形并提取坐标 我正在尝试将该示例改编为不同的形状文件，但出现错误。新的 shapefile 可以从这里下载：https://login.filesanywhere.com/fs/v.aspx?v=8c6d62865f6575ad9ca5

这是我正在使用的新代码：

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)
        )
}

#helper function to rotate the geometry of a spatial dataframe
    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)
    }
    
#load libraries, create play data
    library(sf)
    library(dplyr)
       
    the_lake <- readRDS('vallejo.RDS') ## spatial feature supplied by OP
    the_anchor <- c(5.66e5, 4.21e6) ## lat 38.067690   long -122.244700
    
    
#generate desired spatial features
    the_strips <- 
      get_strips(crs = st_crs(the_lake)) |>
      rotate_feature(angle = 10, anchor = the_anchor)
      the_strips
    
    the_cropped_strips <- the_strips |>  #GETTING AN ERROR HERE....
      st_intersection(the_lake)
      
      Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented,  : 
  Loop 0 is not valid: Edge 7032 has duplicate vertex with edge 7035
In addition: Warning messages:
1: In st_is_longlat(x) :
  bounding box has potentially an invalid value range for longlat data
2: In st_is_longlat(x) :
  bounding box has potentially an invalid value range for longlat data

 ## obtain the corner points (not areas) on intersection:
    the_corners <- 
      st_intersection(
        st_cast(the_strips, 'MULTILINESTRING'),
        st_cast(the_lake, 'MULTILINESTRING'))
         the_corners

Answer 1

您看到的错误是由无效的多边形引起的：

the_lake <- readRDS("vallejo.RDS") 
st_is_valid(the_lake, reason = TRUE)
#> [1] "Loop 0: Edge 7032 has duplicate vertex with edge 7035"

st_make_valid()

会有所帮助，但是上一个问题和当前问题中的形状之间的显着差异是CRS：在投影UTM之前，而现在它是未投影的WGS84。如果您将 cre 设置为 WGS84，当前的

get_strips()

实现不会抱怨，但您还将获得数万十进制的坐标。

因此，您要么必须变换

the_strips

并使用 lon/lat 作为锚点，要么变换当前的形状；对于后者，形状的循环问题自动成为一种非问题（它没有修复，但你将能够摆脱它）：

library(sf)
library(mapview)

# transform to projected CRS, EPSG:26910 (NAD83 / UTM zone 10N)
the_lake <- readRDS("vallejo.RDS") %>% 
  st_transform(26910)

the_anchor <- c(5.66e5, 4.21e6) ## lat 38.067690   long -122.244700

# increase n
the_strips <-
  get_strips(n = 75, crs = st_crs(the_lake)) |>
  rotate_feature(angle = 10, anchor = the_anchor)

the_cropped_strips <- the_strips |>
  st_intersection(the_lake)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

mapview(the_cropped_strips, zcol = "strip_id")

^{创建于 2023-11-19，使用 reprex v2.0.2}