如何使用sf真正计算球形voronoi图？

我想使用世界的球形性质（不是它的投影）制作带有 voronoi 镶嵌的世界地图，类似于使用 D3.js，但使用 R。

据我所知（“再见平坦的地球，欢迎 S2 球面几何”）该sf包现在完全基于该s2包并且应该按照我的需要执行。但我不认为我得到了预期的结果。一个可重现的例子：

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)  

整个南极洲应该比其他两点更靠近墨尔本。我在这里缺少什么？如何计算维诺球体上使用sf？