pacman::p_load(
  concaveman,    # concave hull
  ggalt,         # geom_encircle
  gridExtra,
  INLA,          # inla.nonconvx.hll
  sf,
  spatstat,
  tidyverse
)

Given a set of points, I want to compute a polygon that encompasses the ensemble.

set.seed(20190419)

nclust <- function(x0, y0, radius, n) {
  return(runifdisc(n, radius, centre = c(x0, y0)))
}

pts <- 
  rPoissonCluster(4, 0.2, nclust, radius=0.2, n=20) %>% 
  coords() %>% 
  as.matrix() %>% 
  st_multipoint(dim = "XY") %>% 
  st_sfc() %>% 
  st_cast("POINT") %>% 
  st_sf()
# plot(pts, pch = 20)
pts %>% 
  ggplot() +
  geom_sf()

Convex hull

Simplest solution. However, it cannot identify separate clusters, and can cover vast areas with null density of points.

Note: the convex hull can be computed with grDevices::chull().

However, the package ggalt provides the convenience plotting function geom_encircle() with some nice options like an expansion factor and some slight control on the smoothness of the boundary.

(p_chull <- 
  pts %>%
  ggplot() +
  geom_encircle(
    aes(X, Y),
    data = pts %>% st_coordinates() %>% data.frame(),
    expand = 0.05,
    s_shape = 1.5,
    fill = "gray90"
  ) +
  geom_sf())

Contour of Kernel density estimate

Take a contour plot at some convenient level from a kernel density estimate.

Problem: need to chose a level to make sure all the points are contained. This is costly to do automatically: need to compute several and pick the first that contains all the points.

dens <-
  st_coordinates(pts) %>%
  data.frame() %>%
  with(.,
       MASS::kde2d(
         X,
         Y,
         lims = c(min(X) - 2*sd(X),
                  max(X) + 2*sd(X),
                  min(Y) - 2*sd(Y),
                  max(Y) + 2*sd(Y)),
         h = .5
       )) 

## Manually and carefully chose the level
dens$z <- log(dens$z)
# hist(as.numeric(dens$z))
# plot(st_coordinates(pts), xlim = c(-.5, 1.5), ylim = c(-.5, 1.5))
# contour(dens, add = T)

envelope_ckde <-
  map(
    contourLines(dens, levels = -.1),
    ~ data.frame(.) %>% select(-level) %>% as.matrix %>%
      rbind(., .[1,])  # close polygons
  ) %>% 
  st_polygon() %>% 
  st_sfc(crs = st_crs(pts))
# plot(envelope_ckde)
(p_ckde <- 
  pts %>%
  ggplot() +
  geom_sf(data = envelope_ckde) +
  geom_sf() +
  coord_sf(expand = FALSE)  # remove expansion to panel area
)

Concave hull

The package concaveman ports a very fast 2D concave-hull algorithm from mapbox.

envelope_cm <- 
  concaveman(pts, concavity = 3.5, length_threshold = .1)
(p_cm <- 
  pts %>%
  ggplot() +
  geom_sf(data = envelope_cm) +
  geom_sf() +
  coord_sf(expand = FALSE)  # remove expansion to panel area
)

Here again, an expansion parameter would be desirable in some situations. We can extend a bit beyond the points using a buffer and taking points from the joint polygons as follows.

## This does not extend beyond the points, so I use a buffer and take
## points on the boundaries
envelope_cmb <- 
  pts %>% 
  st_buffer(dist = .05) %>% 
  st_union() %>% 
  st_segmentize(1e3) %>% 
  st_cast("POINT") %>% 
  st_sf() %>% 
  concaveman(
    concavity = 5,
    length_threshold = 1e3
  )
(p_cmb <- 
  pts %>%
  ggplot() +
  geom_sf(data = envelope_cmb) +
  geom_sf() +
  coord_sf(expand = FALSE)  # remove expansion to panel area
)

INLA nonconvex-hull

This provides the most parsimonious and flexible envelope, in my opinion. However, it can’t split the envelope in two disconnected components.

## Using INLA fmesher
envelope_inch <- 
  pts %>% 
  st_coordinates() %>% 
  inla.nonconvex.hull(
    convex = .05,   # buffer distance from points
    concave = .5,   # minimum curvature radius
    resolution = 100
  ) %>% 
  with(., st_polygon(list(rbind(loc, loc[1,])))) %>% 
  st_sfc(crs = st_crs(pts))
(p_inch <- 
  pts %>%
  ggplot() +
  geom_sf(data = envelope_inch) +
  geom_sf() +
  coord_sf(expand = FALSE)  # remove expansion to panel area
)

Side to side comparison