Search code examples
rleafletmapstmapmapview

How to make small polygons from shapefile and extract coordinates


I have a small shapefile here: https://login.filesanywhere.com/fs/v.aspx?v=8c6d62865c626fb4a2ab called bay.RDS

 library(tmap)
 library(leaflet)
 library(mapview)
 bay <- readRDS('bay.RDS')
 mapview(bay)

I am trying to draw 50 meter wide polygons from shore to shore right at the below point and extract coordinates from it. I need multiple polygons next to each other but independently drawn as I need to color code them based on a column of values

pt <- st_sfc(st_point(c(-122.91, 38.17)), crs=4326)
mapview(bay) + mapview(pt, color ='red')

I have included a snapshot below of what I am trying to accomplish. I would appreciate it if someone can show me how I can make just a couple of polygons and extract the coordinates from it so I can merge them to my dataset. Let me know if is not clear

enter image description here


Solution

  • 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.

    1. 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)
            )
    }
    
    1. 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)
        }
    
    1. 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
    
    1. 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')
          )
    
    1. 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))
    
    

    result map

    1. 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