How can I draw a line from a point to a polygon edge and then get the line's length in sf in R?

Question

There are some other posts out there related to this one, such as these: Post 1, Post 2, Post 3. However, none of them deliver what I am hoping for. What I want is to be able to draw a line segment from a specific point (a sampling location) to the edge of a polygon fully surrounding that point (a lake border) in a specific direction ("due south" aka downward). I then want to measure the length of that line segment in between the sampling point and the polygon edge (really, it's only the distance I want, so if we can get the distance without drawing the line segment, so much the better!). Unfortunately, it doesn't seem like functionality to do this already exists within the sf package: See closed issue here.

I suspect, though, that this is possible through a modification of the solution offered here: See copy-pasted code below, modified by me. However, I am pretty lousy with the tools in sf--I got as far as making line segments that just go from the points themselves to the southern extent of the polygon, intersecting the polygon at some point:

library(sf)
library(dplyr)

df = data.frame(
  lon = c(119.4, 119.4, 119.4, 119.5, 119.5),
  lat = c(-5.192,-5.192,-5.167,-5.167,-5.191)
)

polygon <- df %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

plot(polygon)

df2 <- data.frame(lon = c(119.45, 119.49, 119.47),
                  lat = c(-5.172,-5.190,-5.183))

points <- df2 %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("MULTIPOINT")

plot(points, add = TRUE, col = "red")

# Solution via a loop

xmin <- min(df$lat)

m = list()
# Iterate and create lines
for (i in 1:3) {
  m[[i]] = st_linestring(matrix(
    c(df2[i, "lon"],
      df2[i, "lat"],
      df2[i, "lon"],
      xmin),
    nrow = 2,
    byrow = TRUE
  ))
}
test = st_multilinestring(m)

# Result is line MULTILINESTRING object

plot(test, col = "green", add = TRUE)

But now I can't figure out how to use st_intersection or any such function to figure out where the intersection points are. Most of the trouble lies, I think, in the fact that what I'm creating is not an sf object, and I can't figure out how to get it to be one. I assume that, if I could figure out where the segments intersect the polygon (or the most-northern time they do so, ideally), I could somehow measure from the intersection points to the sampling points using a function like st_distance. Since lake polygons are often really complex, though, it's possible a segment will intersect the polygon multiple times (such as if there is a peninsula south of a given point), in which case I figure I can find the "furthest north" intersection point for each sampling point and use that or else take the minimum such distance for each sampling point.

Anyhow, if someone can show me the couple of steps I'm missing, that'd be great! I feel like I'm so close and yet so far...

user_7756889 · Answer

Consider this approach, loosely inspired by my earlier post about lines from points

To make it more reproducible I am using the well known & much loved North Carolina shapefile that ships with {sf} and a data frame of three semi-random NC cities.

What the code does is:

iterates via for cycle over the dataframe of cities
creates a line starting in each city ("observation") and ending on South Pole
intersects the line with dissolved North Carolina
blasts the intersection to individual linestrings
selects the linestring that passes within 1 meter of origin
calculates the lenght via sf::st_lenghth()
saves the the result as a {sf} data frame called res (short for result :)

I have included the actual line in the final object to make the result more clear, but you can choose to omit it.

library(sf)
library(dplyr)
library(ggplot2)

shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%  # included with sf package
  summarise() %>% 
  st_transform(4326) # to align CRS with cities

cities <- data.frame(name = c("Raleigh", "Greensboro", "Plymouth"),
                     x = c(-78.633333, -79.819444, -76.747778),
                     y = c(35.766667, 36.08, 35.859722)) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326)

# a quick overview
ggplot() +
   geom_sf(data = shape) + # polygon of North Carolina
   geom_sf(data = cities, color = "red") # 3 cities

# now here's the action!!!
for (i in seq_along(cities$name)) {
  
  # create a working linestring object
  wrk_line <- st_coordinates(cities[i, ]) %>% 
    rbind(c(0, -90)) %>% 
    st_linestring() %>% 
    st_sfc(crs = 4326) %>% 
    st_intersection(shape) %>% 
    st_cast("LINESTRING") # separate individual segments of multilines
  
  first_segment <- unlist(st_is_within_distance(cities[i, ], wrk_line, dist = 1))  
  
  # a single observation
  line_data <- data.frame(
    name = cities$name[i],
    length = st_length(wrk_line[first_segment]),
    geometry = wrk_line[first_segment]
  )
  
  # bind results rows to a single object
  if (i == 1) {
    res <- line_data
    
  } else {
    res <- dplyr::bind_rows(res, line_data)
    
  } # /if - saving results
  
  
} # /for

# finalize results
res <- sf::st_as_sf(res, crs = 4326) 
  
# result object    
res
# Simple feature collection with 3 features and 2 fields
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: -79.81944 ymin: 33.92945 xmax: -76.74778 ymax: 36.08
# Geodetic CRS:  WGS 84
#         name        length                       geometry
# 1    Raleigh 204289.21 [m] LINESTRING (-78.63333 35.76...
# 2 Greensboro 141552.67 [m] LINESTRING (-79.81944 36.08...
# 3   Plymouth  48114.32 [m] LINESTRING (-76.74778 35.85...

# a quick overview of the lines
ggplot() +
  geom_sf(data = shape) + # polygon of North Carolina
  geom_sf(data = res, color = "red") # 3 lines

How can I draw a line from a point to a polygon edge and then get the line's length in sf in R?

1 Answers