Does the order of locations change when making a voronoi map with leaft and sf?

I am making a voronoi map in R using the packages leaflet and sf as follows:

library(leaflet)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

long   <- c(4.35556 ,  5.83745,  4.63683 ,  6.06389,  6.41111,  5.639722)
lat    <- c(52.00667, 53.09456, 52.38084 , 52.475  , 52.15917, 53.440278)
labs   <- c("Delft" , "Grouw" , "Haarlem", "Hattem", "Lochem", "Hollum" )
colors <- c("red"   , "orange", "yellow" , "green" , "blue"  , "purple" )

df <- data.frame(ID = labs, X = long, Y = lat)
points <- st_geometry(st_as_sf(df, coords = c("X", "Y")))
points <- st_union(points)

NL <- ne_countries(country = c('netherlands'), scale = 'medium', returnclass = 'sf')

polys <- points       %>% 
  st_voronoi()        %>%
  st_cast()           %>%
  st_set_crs(., 4326) %>%
  st_intersection(y = NL)

leaflet() %>%
addProviderTiles(providers$OpenStreetMap.Mapnik) %>%
addPolygons     (data        = polys,
                 fillColor   = colors,
                 fillOpacity = 1,
                 weight      = 0.5, 
                 color       = "black") %>%
addCircleMarkers(lng         = long,
                 lat         = lat,
                 label       = labs, 
                 color       = "black",
                 radius      = 5,
                 weight      = 1,
                 fill        = TRUE,
                 fillColor   = colors,
                 fillOpacity = 1)

In the resulting map the colors of the dots are correct, but the colors of the polygons are not correct. I guess something has changed in order of the locations in 'polys', but I am puzzled about this. Any suggestions how to solve this?

Solution

st_voronoi() indeed appears not to keep the order of input points in the resulting polygons. You may use st_intersects() to find out which polygon belongs to which point and reorder polys accordingly.

First store a copy of points before applying st_union() and set them the same CRS as polys will have, so thatst_intersects() works later on. I.e., insert this before the points <- st_union(points) line:

points0 <- st_set_crs(points, 4326)

Then, after creating polys, reorder them like this:

polys <- polys[unlist(st_intersects(points0, polys))]

If some point is located outside the area of Netherlands (as provided by ne_countries()) the matching of points to polygons has to be done before intersecting of polys and NL. So in the original code the polys <- points... will be replaced with:

polys <- points       %>% 
  st_voronoi()        %>%
  st_cast()           %>%
  st_set_crs(., 4326) 
polys <- polys[unlist(st_intersects(points0, polys))]
polys <- st_intersection(polys, NL)