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