How to truly calculate a spherical voronoi diagram using sf?

Question

I want to make a world map with a voronoi tessellation using the spherical nature of the world (not a projection of it), similar to this using D3.js, but with R.

As I understand ("Goodbye flat Earth, welcome S2 spherical geometry") the sf package is now fully based on the s2 package and should perform as I needed. But I don't think that I am getting the results as expected. A reproducible example:

library(tidyverse)
library(sf)
library(rnaturalearth)
library(tidygeocoder)

# just to be sure
sf::sf_use_s2(TRUE)

# download map 
world_map <- rnaturalearth::ne_countries(
                scale = 'small', 
                type = 'map_units',
                returnclass = 'sf')

# addresses that you want to find lat long and to become centroids of the voronoi tessellation 
addresses <- tribble(
~addr,
"Juneau, Alaska" ,
"Saint Petersburg, Russia" ,
"Melbourne, Australia" 
)

# retrive lat long using tidygeocoder
points <- addresses %>% 
          tidygeocoder::geocode(addr, method = 'osm')

# Transform lat long in a single geometry point and join with sf-base of the world
points <- points %>% 
          dplyr::rowwise() %>% 
          dplyr::mutate(point = list(sf::st_point(c(long, lat)))) %>% 
          sf::st_as_sf() %>% 
          sf::st_set_crs(4326)

# voronoi tessellation
voronoi <- sf::st_voronoi(sf::st_union( points ) ) %>% 
     sf::st_as_sf() %>% 
     sf::st_set_crs(4326)

# plot
ggplot2::ggplot() +
    geom_sf(data = world_map,
            mapping = aes(geometry = geometry), 
            fill = "gray95") +
    geom_sf(data = points,
            mapping = aes(geometry = point),
            colour = "red") +
    geom_sf(data = voronoi,
            mapping = aes(geometry = x),
            colour = "red",
            alpha = 0.5)  

The whole Antarctica should be closer to Melbourne than to the other two points. What am I missing here? How to calculate a voronoi on a sphere using sf?

Answer 1

Here's a method that builds on Stéphane Laurent's approach, but outputs sf objects.

Let us obtain an sf object of all the world capitals:

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

capitals <- do.call(rbind,
  subset(maps::world.cities, capital == 1, select = c("long", "lat")) |>
  as.matrix() |>
  asplit(1) |>
  lapply(st_point) |>
  lapply(st_sfc) |>
  lapply(st_sf, crs = 'WGS84')) |>
  `st_geometry<-`('geometry') |>
  cbind(city = subset(maps::world.cities, capital == 1, select = c("name")))

capitals
#> Simple feature collection with 230 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -176.13 ymin: -51.7 xmax: 179.2 ymax: 78.21
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>           name               geometry
#> 1       'Amman    POINT (35.93 31.95)
#> 2    Abu Dhabi    POINT (54.37 24.48)
#> 3        Abuja      POINT (7.17 9.18)
#> 4        Accra      POINT (-0.2 5.56)
#> 5    Adamstown  POINT (-130.1 -25.05)
#> 6  Addis Abeba     POINT (38.74 9.03)
#> 7        Agana   POINT (144.75 13.47)
#> 8      Algiers     POINT (3.04 36.77)
#> 9        Alofi POINT (-169.92 -19.05)
#> 10   Amsterdam     POINT (4.89 52.37)

And our world map:

world_map <- rnaturalearth::ne_countries(
  scale = 'small', 
  type = 'map_units',
  returnclass = 'sf')

Now we use Stéphane Laurent's approach to tesselate the sphere, but then reverse the projection back into spherical co-ordinates. This allows translation back to sf, though we have to be careful to split any objects that "wrap around" the 180/-180 longitude line:

voronoi <- capitals %>%
  st_coordinates() %>%
  `*`(pi/180) %>%
  cbind(1) %>%
  pracma::sph2cart() %>%
  sphereTessellation::VoronoiOnSphere() %>%
  lapply(\(x) rbind(t(x$cell), t(x$cell)[1,])) %>%
  lapply(\(x) {
    n <- nrow(x) - 1
    lapply(seq(n), function(i) {
      a <- approx(x[i + 0:1, 1], x[i + 0:1, 2], n = 1000)
      b <- approx(x[i + 0:1, 1], x[i + 0:1, 3], n = 1000)
      d <- cbind(a$x, a$y, b$y) |> pracma::cart2sph() 
      d <- d[,1:2] * 180/pi
      if(max(abs(diff(d[,1]))) > 180) {
        s <- which.max(abs(diff(d[,1])))
        d <- list(d[1:s, ], d[(s+1):nrow(d),])
      }
      d })}) |>
  lapply(\(x) {
    st_geometrycollection(lapply(x, \(y) {
    if(class(y)[1] == "list") {
      st_multilinestring(y)
      } else {
        st_linestring(y)
      }}))}) %>%
  lapply(st_sfc) %>%
  lapply(st_sf, crs = 'WGS84') %>%
  {do.call(rbind, .)} %>%
  `st_geometry<-`('geometry')

Now we have our Voronoi grid as an sf object, so we can plot it using ggplot:

library(ggplot2)

ggplot() +
  geom_sf(data = world_map, fill = "cornsilk", color = NA) +
  geom_sf(data = voronoi, color = "gray40") +
  geom_sf(data = capitals, color = "black", size = 0.2) + 
  coord_sf(crs = "ESRI:53011") +
  theme(panel.background = element_rect(fill = "lightblue"))

---

Addendum

Although the above solution works for drawing tesselations over the whole globe, if we want to get polygons of land areas only, we can do it as follows:

First, we make a union of all land masses from our world map

wm <- st_make_valid(world_map) |> st_union()

Now we get the co-ordinates of the vertices of our Voronoi tiles:

pieces <- capitals %>%
  st_coordinates() %>%
  `*`(pi/180) %>%
  cbind(1) %>%
  pracma::sph2cart() %>%
  sphereTessellation::VoronoiOnSphere() %>%
  lapply(\(x) rbind(t(x$cell), t(x$cell)[1,])) %>%
  lapply(pracma::cart2sph) %>%
  lapply(\(x) x[,1:2] * 180/pi)

Now we need to find tiles that span the -180 / 180 line:

complete <- pieces %>% sapply(\(x) abs(diff(c(min(x[,1]), max(x[,1])))) < 180)

We now split these and turn them into multipolygons, finding their intersection with the world map:

orphans <- pieces[!complete] %>%
  lapply(\(x) {x[,1] + 180 -> x[,1]; x}) %>%
  lapply(\(x) st_polygon(list(x)) |> st_sfc(crs = "WGS84")) %>%
  lapply(\(x) {
    west <- st_intersection(x, matrix(c(-180, -0.001, -0.001, -180, -180, 
                                 -89, -89, 89, 89, -89), ncol = 2) |>
                                list() |> st_polygon() |> st_sfc(crs = "WGS84"))
    east <- st_intersection(x, matrix(c(0, 180, 180, 0, 0, 
                                        -89, -89, 89, 89, -89), ncol = 2) |>
                              list() |> st_polygon() |> st_sfc(crs = "WGS84"))
    west <- st_coordinates(west)[,1:2]
    east <- st_coordinates(east)[,1:2]
    west[,1] <- west[,1] + 180
    east[,1] <- east[,1] - 180
    w <- st_polygon(list(west)) |> st_sfc(crs = "WGS84") |> st_intersection(wm)
    e <- st_polygon(list(east)) |> st_sfc(crs = "WGS84") |> st_intersection(wm)
    st_combine(st_union(e, w))
  }) %>%
  lapply(st_sf) %>%
  lapply(\(x) { if(nrow(x) > 0) { st_segmentize(x, 100000) } else {
    st_point(matrix(c(0, 0), ncol = 2)) |> 
      st_sfc(crs = "WGS84") |> st_sf()
  }
  }) %>% 
  lapply(\(x) `st_geometry<-`(x, 'geometry')) %>%
  {do.call(rbind, .)} %>%
  cbind(city = capitals$name[!complete])

We can do intersections for the non-wraparound tiles like this:

non_orphans <- pieces %>%
  subset(complete) %>%
  lapply(list) %>%
  lapply(st_polygon) %>%
  lapply(st_sfc, crs = "WGS84") %>%
  lapply(st_intersection, y = wm) %>%
  lapply(st_sf) %>%
  lapply(\(x) { if(nrow(x) > 0) { st_segmentize(x, 100000) } else {
    st_point(matrix(c(0, 0), ncol = 2)) |> 
      st_sfc(crs = "WGS84") |> st_sf()
  }
  }) %>%
  lapply(\(x) `st_geometry<-`(x, 'geometry')) %>%
  {do.call(rbind, .)} %>%
  cbind(city = capitals$name[complete])

Finally, we bind all these together into a single sf object:

voronoi <- rbind(orphans, non_orphans)
voronoi <- voronoi[!st_is_empty(voronoi),]
voronoi <- voronoi[sapply(voronoi$geometry, \(x) class(x)[2] != "POINT"),]

Now we're ready to plot. Let's define a palette function that gives results similar to your linked example:

f <- colorRampPalette(c("#dae7b4", "#c5b597", "#f3dca8", "#b4b6e7", "#d6a3a4"))

We'll also create a background "globe" and a smoothed grid to draw over our map, as in the example:

grid <- lapply(seq(-170, 170, 10), \(x) rbind(c(x, -89), c(x, 0), c(x, 89))) |>
  lapply(st_linestring) |>
  lapply(\(x) st_sfc(x, crs = "WGS84")) |>
  lapply(\(x) st_segmentize(x, dfMaxLength = 100000)) |>
  c(
    lapply(seq(-80, 80, 10), \(x) rbind(c(-179, x), c(0, x), c(179, x))) |>
      lapply(st_linestring) |>
      lapply(\(x) st_sfc(x, crs = "WGS84"))
  ) |>
  lapply(st_sf) |>
  lapply(\(x) `st_geometry<-`(x, 'geometry')) %>%
  {do.call(rbind, .)}

globe <- st_polygon(list(cbind(c(-179, 179, 179, -179, -179), 
                               c(-89, -89, 89, 89, -89)))) |> 
  st_sfc(crs = "WGS84") |> 
  st_segmentize(100000)

The final result is a faithful sf version of the linked example:

ggplot() +
  geom_sf(data = globe, fill = "#4682b4", color = "black") +
  geom_sf(data = voronoi, color = "black", aes(fill = city)) +
  geom_sf(data = capitals, color = "black", size = 1) + 
  geom_sf(data = grid, color = "black", linewidth = 0.2) +
  coord_sf(crs = "ESRI:53011") +
  scale_fill_manual(values = f(nrow(voronoi))) +
  theme(panel.background = element_blank(),
        legend.position = "none",
        panel.grid = element_blank())

^{Created on 2023-06-24 with reprex v2.0.2}

Answer 2

With sf I don't know, but you can use the sphereTessellation package instead.

library(pracma) # for the sph2cart function
library(maps)
data(world.cities)
data(worldMapEnv)

world <- map("world", plot = FALSE)
countries <- sph2cart(cbind(world$x*pi/180, world$y*pi/180, 1))

capitals_ll <- as.matrix(
    subset(world.cities, capital == 1, select = c("long", "lat"))
) * pi / 180
capitals <- sph2cart(cbind(capitals_ll, 1))


library(rgl)
open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.7)
sphereMesh <- Rvcg::vcgSphere()
shade3d(sphereMesh, color = "cyan", polygon_offset = 1)
lines3d(countries)
points3d(capitals, size = 4, color = "red")
snapshot3d("world.png", webshot = FALSE)

library(sphereTessellation)
vor <- VoronoiOnSphere(capitals)
open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.7)
plotVoronoiOnSphere(
    vor, colors = "white", edges = TRUE, ecolor = "green",
    sites = TRUE, scolor = "red", polygon_offset = 1
)
lines3d(countries)
snapshot3d("world_voronoi.png", webshot = FALSE)

But I don't know how we could clip the Voronoï cells to the countries.

Answer 3

(This answer doesn't tell you how to do it, but does tell you what's going wrong.)

When I ran this code I got

Warning message: In st_voronoi.sfc(sf::st_union(points)) : st_voronoi does not correctly triangulate longitude/latitude data

From digging into the code it looks like this is a known limitation. Looking at the C++ code for CPL_geos_voronoi, it looks like it directly calls a GEOS method for building Voronoi diagrams. It might be worth opening an sf issue to indicate that this is a feature you would value (if no-one tells the developer that particular features would be useful, they don't get prioritized ...) It doesn't surprise me that GEOS doesn't automatically do computations that account for spherical geometry. Although the S2 code base mentions Voronoi diagrams in a variety of places, it doesn't look like there is a drop-in replacement for the GEOS algorithm ... there are a variety of implementations in other languages for spherical Voronoi diagrams (e.g. Python), but someone would probably have to port them to R (or C++) ...

If I really needed to do this I would probably try to figure out how to call the Python code from within R (exporting the data from sf format to whatever Python needs, then re-importing the results into an appropriate sf format ...)

Printing the code for sf:::st_voronoi.sfc:

function (x, envelope = st_polygon(), dTolerance = 0, bOnlyEdges = FALSE) 
{
    if (compareVersion(CPL_geos_version(), "3.5.0") > -1) {
        if (isTRUE(st_is_longlat(x))) 
            warning("st_voronoi does not correctly triangulate longitude/latitude data")
        st_sfc(CPL_geos_voronoi(x, st_sfc(envelope), dTolerance = dTolerance, 
            bOnlyEdges = as.integer(bOnlyEdges)))
    }
    else stop("for voronoi, GEOS version 3.5.0 or higher is required")
}

In other words, if the GEOS version is less than 3.5.0, the operation fails completely. If it is >= 3.5.0 (sf:::CPL_geos_version() reports that I have version 3.8.1), and long-lat data is being used, a warning is supposed to be issued (but the computation is done anyway).

The first time I ran this I didn't get the warning; I checked and options("warn") was set to -1 (suppressing warnings). I'm not sure why — running from a clean session did give me the warning. Maybe something in the pipeline (e.g. rnaturalearth telling me I needed to install the rnaturalearthdata package) accidentally set the option?