# Distance from multilinestring with sf

I have a `MULTILINESTRING` `sf` object, corresponding to rivers and I would like to create a `raster` containing the distance values to the closest river.

I managed to get this with a single river but when I have several of them, I do not obtain a convincing result. I provide hereafter a reproducible example with the Seine river to illustrate my problem.

Example with one segment only

Here, I use only one department as a region of interest and crop the seine river to this department. The following code is inspired by https://dominicroye.github.io/en/2019/calculating-the-distance-to-the-sea-in-r/. It basically draws a grid having the same extent than the zone of interest and then computes for each cell of this grid the distance to the closest river.

``````# Import Seine river
library(spData)
data("seine")

# Import one department
library(maps)
map("france", namesonly = TRUE, plot = FALSE) # get names of French departments
sm <- map("france", regions = "Seine-Maritime", fill = FALSE, col = "black",
plot = FALSE)

# Convert to sf objects
library(sf)
sm_sf <- st_as_sf(sm)
seine_sf <- st_as_sf(seine)

# Crop seine to Seine-Maritime department
seine_int <- seine_sf %>%
st_transform(crs = st_crs(sm_sf)) %>% # same CRS
st_intersection(sm_sf) # crop

# Draw grid
sm_grid <- st_make_grid(sm_sf, cellsize = 0.05, what = "centers")
# Crop the grid to Seine-Maritime department
sm_grid_int <- st_intersection(sm_grid, sm_sf)
# Compute distance to Seine river
dist_seine <- st_distance(seine_int, sm_grid_int)
# Conversion to dataframe
dist_seine_df <- data.frame(dist_seine = as.vector(dist_seine)/1000,
st_coordinates(sm_grid_int))

# Plot
library(ggplot2)
library(cowplot)
library(RColorBrewer)
plot_grid(
ggplot(sm_sf) +
geom_sf(fill = "white") +
geom_sf(data = seine, color = "dodgerblue"),
ggplot(sm_sf) +
geom_sf(fill = "white") +
geom_sf(data = seine_int, color = "dodgerblue"),
ggplot(sm_grid_int) +
geom_sf(),
ggplot(dist_seine_df, aes(X, Y, fill = dist_seine)) +
geom_tile() +
geom_sf(data = sm_sf, inherit.aes = FALSE, fill = NA, size = 1) +
geom_sf(data = seine_int, inherit.aes = FALSE, color = "white") +
labs(fill = "Distance to water (km)") +
theme_void() +
theme(legend.position = "bottom"),
nrow = 2, rel_widths = rep(1, 4), rel_heights = rep(1, 4), align = "vh")
`````` The obtained raster seems legit.

Example with several segments

Now I use the same code as above except that I do not limit it to one department but use all French departments instead.

``````fr_departments <- map("france", fill = FALSE, col = "black", plot = FALSE)

fr_departments_sf <- st_as_sf(fr_departments)
seine_sf <- st_as_sf(seine)

seine_int <- seine_sf %>%
st_transform(crs = st_crs(fr_departments_sf)) %>% # same CRS
st_intersection(fr_departments_sf) # crop

fr_departments_grid <- st_make_grid(fr_departments_sf, cellsize = 0.5,
what = "centers")
# fr_departments_grid_int <- st_intersection(fr_departments_grid,
#                                            fr_departments_sf)
# https://github.com/r-spatial/sf/issues/347
fr_departments_grid_int <- st_intersection(fr_departments_grid,
st_buffer(fr_departments_sf, 0))

# Compute distance to Seine river
dist_seine <- st_distance(seine_int, fr_departments_grid_int)
# RDS object saved: cellsize = 0.1 for the grid
dist_seine_df <- data.frame(dist_seine = as.vector(dist_seine)/1000,
st_coordinates(fr_departments_grid_int))

# Plot
plot_grid(
ggplot(fr_departments_sf) +
geom_sf(fill = "white") +
geom_sf(data = seine_int, color = "dodgerblue"),
NULL,
ggplot(fr_departments_grid_int) +
geom_sf(),
ggplot(dist_seine_df, aes(X, Y, fill = dist_seine)) +
geom_tile() +
geom_sf(data = fr_departments_sf, inherit.aes = FALSE, fill = NA, size = 1) +
geom_sf(data = seine_int, inherit.aes = FALSE, color = "white") +
labs(fill = "Distance to water (km)") +
theme_void() +
theme(legend.position = "bottom"),
nrow = 2, rel_widths = rep(1, 4), rel_heights = rep(1, 4), align = "vh")
`````` Here we can see that the obtained raster gives weird output with remote cells having short distance from the river.

What did I do wrong with this code?

Side question

`st_intersection` can be very slow to run depending on the resolution of the grid. Is there a faster `sf` solution to this?

The answer was actually straightforward. With a `MULTILINESTRING` object, a distance is computed between each raster cell of the grid and each segment. If interested into the closest distance, the minimum distance has to be extracted.

Therefore:

``````# Compute distance to Seine river
dist_seine <- st_distance(seine_int, fr_departments_grid_int)
# Only minimum distance
dist_seine <- apply(dist_seine, 2, min)
``````

Edit

`union` the features before computing the distance avoids computing the minimum distance and gives similar results. However, `microbenchmark` suggests that the computation time is slightly slower in that case.

``````library(microbenchmark)
comp <- microbenchmark(
apply_min = {dist_seine <- st_distance(seine_int, fr_departments_grid_int)
dist_seine <- apply(dist_seine, 2, min)},

union = dist_union_seine <- st_distance(st_union(st_geometry(seine_int)),
fr_departments_grid_int))
library(ggplot2)
autoplot(comp)
`````` • You can `union` your features to get one feature, and then `st_distance` returns single-row matrix, eg: `dist_union_seine <- st_distance(st_union(st_geometry(seine_int)), fr_departments_grid_int)` and there's no need to sweep the minimum. Might be quicker... Nov 29 '19 at 23:14
• Thank you for the suggestion, `union` the features does give similar results but it appears to be slightly slower. I edited the answer so both options are visible. Dec 2 '19 at 10:04
• Thanks for considering and testing! Might be quicker if you have many many segments and many many points - should also take less memory, since it doesn't have to create an NxM ful distance matrix. Dec 2 '19 at 10:06
• I was also wondering, that's why I put both options in the answer. Dec 2 '19 at 10:07