3

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") +
    scale_fill_gradientn(colours = rev(brewer.pal(9, "Blues"))) +
    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")

enter image description here

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 <- readRDS("D:/PIERRE_DENELLE/Stackoverflow/world_dist_rivers.rds")
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") +
    scale_fill_gradientn(colours = rev(brewer.pal(9, "Blues"))) +
    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")

enter image description here

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?

2

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)

which gives: enter image description here

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)

enter image description here

| improve this answer | |
  • 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... – Spacedman 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. – P. Denelle 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. – Spacedman Dec 2 '19 at 10:06
  • I was also wondering, that's why I put both options in the answer. – P. Denelle Dec 2 '19 at 10:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.