How to set an anchor and create polygons from shapefile

Question

I am following up on the example here: How to make small polygons from shapefile and extract coordinates I am trying to adapt that example to a different shapefile but I am getting errors. The new shapefile can be downloaded from here: https://login.filesanywhere.com/fs/v.aspx?v=8c6d62865f6575ad9ca5

This is the new code that I am using:

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 error you see is caused by the input data, polygon is invalid:

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() would help, BUT the significant difference between shapes in your previous question and the current one is the CRS: before it was projected UTM while now it's unprojected WGS84. The current get_strips() implementation will not complain if you set crs to WGS84, but you'll also have coordinates in the rage of tens of thousands decimal degrees; that's what those bounding box has potentially an invalid value range for longlat data warnings are about.

So you either have to transform the_strips and use lon/lat for the anchor or transform your current shape; with the latter, the loop-issue of the shape automgically becomes kind of a non-issue (it's not fixed but you'll be able to get away with it):

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