One approach, using {sf} and {dplyr}:
edit this solution was substantially rewritten following the helpful comments by @Chris and @Salvador.
TLDR: below solution creates a horizontal baseline, buffers this baseline to result in a horizontal strip of chosen width, creates adjacent clones of this strip, combines them into a multipolygon, and rotates and shifts the mulitpolygon to expand southwards from a chosen anchor.
- helper function
get_strips(...) to return a multipolygon of adjacent 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)
)
}
- 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(file.choose()) ## spatial feature supplied by OP
the_anchor <- c(5.08e5, 4.224e6) ## anchor strip grid here
- generate desired spatial features
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')
)
- inspect
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))
- extract corners per polygon:
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
