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

(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?