Can I bound an st_distance call by a polygon?

Question

I have seen similar posts on this topic (see, for example, here and here) but not one that is specific to the sf-tidyverse ecosystem.

I have a series of lakes, a series of sample points within each lake, and a "focal point" in each lake that represents where a boat launch is. I want to calculate the "boater's shortest travel distance" to each sample point from the boat launch. However, I want to somehow "bound" these distances using the lake polygon such that distances cannot be calculated across land. I could imagine this being done by having the "straight line" snake along the lake edge for as long as needed before it can resume being a straight line. I could also imagine the "straight line" being decomposed into line segments that bend around any intervening land. I'm open to a variety of implementations!

I have seen elsewhere (such as here) the idea of converting the bounding polygon to a raster, making the water one value and the land another, much higher value, and then doing a "least-cost distance," where the cost of going over land is prohibitive. However, I don't know how one would actually do this in the raster/sf ecosystem.

Here's the code I used to make this map:

library(dplyr)
library(sf)
library(ggplot2)
Moose.ssw = sswMN.sf %>% filter(lake == "Moose")
Moose.lake = MN_lakes4 %>% filter(str_detect(map_label, "Moose")) %>% filter(cty_name == "Beltrami") 
Moose.access = try3 %>% filter(LAKE_NAME == "Moose") %>%  filter(COUNTYNAME == "Beltrami")
Moose.box = st_bbox(Moose.ssw)

ggplot() +
  geom_sf(data = Moose.lake, color="lightblue") + 
  geom_sf(data = Moose.access, color = "red", size = 2) +
  geom_sf(data = Moose.ssw, mapping = aes(color= Nitellopsis_obtusa_n)) +
  coord_sf(xlim = c(Moose.box[1], Moose.box[3]), ylim = c(Moose.box[2], Moose.box[4]))

And here's some code that calculates straight-line distances splendidly--maybe it can be wrappered somehow?

Moose.pt.dists = st_distance(Moose.access, Moose.ssw, by_element = TRUE)

Files needed to make the data objects referenced above can be pulled from my Github page (they are files produced via dput(). Link to my Github.

I am a solid R programmer but I am new to geospatial work, so if I could even just be pointed in a fruitful direction, I may be able to program my own way out of this!

Answer 1

Here's a solution using sfnetworks, which fits in with the tidyverse well.

The code below should

• regularly sample the bounding box of the lake (creates evenly-spaced points)
• connect them to their closest neighbors
• remove the connections that cross land
• find the shortest path from the boat launch to the sample location(s) by following the lines that remain.

The results are not exact, but are pretty close. You can increase precision by increasing the number of sampled points. 1000 points are used below.

library(tidyverse)
library(sf)
library(sfnetworks)
library(nngeo)

# set seed to make the script reproducible,
#  since there is random sampling used
set.seed(813)

# Getting your data:
x <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
moose_lake <- x[5,]
boat_ramp <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
sample_locations <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
sample_bbox <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")


# sample the bounding box with regular square points, then connect each point to the closest 9 points
#  8 should've worked, but left some diagonals out.
sq_grid_sample <- st_sample(st_as_sfc(st_bbox(moose_lake)), size = 1000, type = 'regular') %>% st_as_sf() %>%
  st_connect(.,.,k = 9)

# remove connections that are not within the lake polygon
sq_grid_cropped <- sq_grid_sample[st_within(sq_grid_sample, moose_lake, sparse = F),]

# make an sfnetwork of the cropped grid
lake_network <- sq_grid_cropped %>% as_sfnetwork()

# find the (approximate) distance from boat ramp to point 170 (far north)
pt170 <- st_network_paths(lake_network, 
                          from = boat_ramp,
                          to = sample_locations[170,]) %>%
  pull(edge_paths) %>%
  unlist()

lake_network %>% 
  activate(edges) %>%
  slice(pt170) %>%
  st_as_sf() %>%
  st_combine() %>%
  st_length()
#> 2186.394 [m]

Looks like it is about 2186m from the boat launch to sample location 170 at the far north end of the lake.

# Plotting all of the above, path from boat ramp to point 170:

ggplot() + 
  geom_sf(data = sq_grid_sample, alpha = .05) +
  geom_sf(data = sq_grid_cropped, color = 'dodgerblue') + 
  geom_sf(data = moose_lake, color = 'blue', fill = NA) +
  geom_sf(data = boat_ramp, color = 'springgreen', size = 4) + 
  geom_sf(data = sample_locations[170,], color = 'deeppink1', size = 4) +
  geom_sf(data = lake_network %>% 
            activate(edges) %>%
            slice(pt170) %>%
            st_as_sf(),
          color = 'turquoise',
          size = 2) +
  theme_void()

Though only the distance from the boat launch to one sample point is found and plotted above, the network is there to find all of the distances. You may need to use *apply or purrr, and maybe a small custom function to find the 'one-to-many' distances of the launch to all sample points.

This page on sfnetworks will be helpful in writing the one-to-many solution.

edit:

To find all distances from the boat launch to the sample points:

st_network_cost(lake_network,
                from = boat_ramp,
                to = sample_locations)

This is much faster for me than a for loop or the sp solution posted below. Some code in the sampling may need to be adjusted since the st_network_cost will remove any duplicate distances. The sample_locations (or Moose.ssw in @bajcz answer) will need to be cleaned to remove duplicate points as well. There are 180 unique points of the 322 rows.