Buffer areas for nearest neighbours or national marine areas delineation

A few weeks ago, Christina Buelow asked on Twitter how she could create the polygons of coastal waters using the land polygons, and making sure there is no overlap in buffer areas. This would make it possible to know which land area is closest to any place in the sea. Let’s explore what solution emerged from the discussions.

Create buffer area around adjacent polygons

As always when trying to figure out a solution to a specific problem, let’s start with a reproducible example, commonly called “reprex” among R users. This means, we create the minimal code that allows anyone to execute the code, face the same problem and show everything we tested. In the present case, I couldn’t find enough time to reach a solution during the discussions, so I prepared a gist showing what I tried, what point I reached, and more importantly what I plan to achieve.

First, we load some packages

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

my_blue <- "#1e73be"

Then we download some data

# Define where to save the dataset
extraWD <- "."
# Get some data available to anyone
if (!file.exists(file.path(extraWD, "departement.zip"))) {
  githubURL <- "https://github.com/statnmap/blog_tips/raw/master/2018-07-14-introduction-to-mapping-with-sf-and-co/data/departement.zip"
  download.file(githubURL, file.path(extraWD, "departement.zip"))
  unzip(file.path(extraWD, "departement.zip"), exdir = extraWD)
}

I reduce the dataset to better see what I am doing. Here, I filter the Brittany area.
As you can see, there are 4 departments in this region, thus 4 spatial features.

departements_l93 <- read_sf(dsn = extraWD, layer = "DEPARTEMENT")

# Reduce dataset
bretagne_l93 <- departements_l93 %>%
  filter(NOM_REG == "BRETAGNE")

bretagne_l93
## Simple feature collection with 4 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 99217.1 ymin: 6704195 xmax: 400572.3 ymax: 6881964
## projected CRS:  RGF93_Lambert_93
## # A tibble: 4 x 12
##   ID_GEOFLA CODE_DEPT NOM_DEPT CODE_CHF NOM_CHF X_CHF_LIEU Y_CHF_LIEU X_CENTROID
## * <chr>     <chr>     <chr>    <chr>    <chr>        <int>      <int>      <int>
## 1 DEPARTEM… 35        ILLE-ET… 238      RENNES      351870    6789646     326237
## 2 DEPARTEM… 22        COTES-D… 278      SAINT-…     274947    6839206     260216
## 3 DEPARTEM… 29        FINISTE… 232      QUIMPER     171326    6789950     115669
## 4 DEPARTEM… 56        MORBIHAN 260      VANNES      267889    6744074     256790
## # … with 4 more variables: Y_CENTROID <int>, CODE_REG <chr>, NOM_REG <chr>,
## #   geometry <MULTIPOLYGON [m]>

Now, I can plot the output of buffer areas around separate polygons.

# Create buffer
bretagne_buffer_l93 <- bretagne_l93 %>% 
  st_buffer(
  dist = 
    units::set_units(10, km)
) 

ggplot() +
  geom_sf(data = bretagne_buffer_l93, fill = my_blue, alpha = 0.2, colour = my_blue) +
  geom_sf(data = bretagne_l93, fill = NA) +
  theme_bw()

You can see there are overlapping areas.
This is where Voronoï polygons come to help us. From here, I will take some ideas that Christina gathered during her Twitter discussions. Her complete code for coastal buffers is here.

# Create a merged region entity
bretagne_union_l93 <- bretagne_l93 %>% 
  summarise()

ggplot() + 
  geom_sf(data = bretagne_union_l93, colour = my_blue, fill = NA)

# Create a doughnut for regional seas areas, 30km off coasts
bretagne_seas_l93 <- bretagne_union_l93 %>% 
  st_buffer(
    dist = units::set_units(30, km)
  ) %>% 
  st_cast()

ggplot() + 
  geom_sf(data = bretagne_seas_l93, colour = my_blue, fill = NA)

# Remove inside terrestrial parts
bretagne_donut_l93 <- bretagne_seas_l93 %>% 
  st_difference(bretagne_union_l93) %>% 
  st_cast()

ggplot() + 
  geom_sf(data = bretagne_donut_l93, colour = my_blue, fill = my_blue, alpha = 0.2)

Let’s divide this area in multiples polygons

We use the coastlines as a basis.

# density <- units::set_units(0.1, 1/km)

# First merge everything and transform as lines
bretagne_points_l93 <- bretagne_union_l93 %>% 
  # transform polygons as multi-lines
  st_cast("MULTILINESTRING") %>% 
  # transform multi-lines as separate unique lines
  st_cast("LINESTRING") %>% 
  # transform as points, 0.1 per km (= 1 every 10 km)
  # Choose density according to precision needed
  st_line_sample(density = units::set_units(0.1, 1/km)) %>% 
  st_sf() #%>% 
  # remove empty geometry (small islands where sample is not possible with this density)
  # filter(!st_is_empty(geometry))
  
ggplot() + 
  geom_sf(data = bretagne_points_l93, colour = my_blue, fill = NA)

# Create voronoi polygons around points
bretagne_voronoi_l93 <- bretagne_points_l93 %>% 
  st_combine() %>% 
  st_voronoi() %>% 
  st_cast()

ggplot() + 
  geom_sf(data = bretagne_voronoi_l93, colour = my_blue, fill = NA)

Join voronoi and donut polygon

bretagne_voronoi_donut_l93 <- bretagne_voronoi_l93 %>% 
  st_intersection(bretagne_donut_l93) %>% 
  st_cast() %>% 
  st_sf()

ggplot() + 
  geom_sf(data = bretagne_voronoi_donut_l93, colour = my_blue, fill = NA)

Join voronoi-donut polygons to original dataset to intersect with regions

We realise sort of a left_join based on spatial attributes using st_join(). This will keep the original geometry of our voronoi-donut polygons, while retrieving variables of joined terrestrial dataset when there is intersection.

# Equivalent to left_join
bretagne_vd_region_l93 <- bretagne_voronoi_donut_l93 %>% 
  st_join(bretagne_l93)

ggplot() + 
  geom_sf(data = bretagne_vd_region_l93,
          aes(colour = NOM_DEPT, fill = NOM_DEPT),
          alpha = 0.25)

# Interactive view
mapview::mapview(bretagne_vd_region_l93, zcol = "NOM_DEPT")