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

Question

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?

Robert Hacken · Accepted Answer · 2022-09-20T17:53:35.543

1

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)

edited Sep 20 '22 at 17:53

answered Sep 15 '22 at 16:21

Robert Hacken

3,878
1
13
15

2

The help file for `st_voronoi` has this in its examples and has saved me several times.... – Spacedman Sep 15 '22 at 17:06
2

@Spacedman Oh, indeed! And switching polygons and points in `st_intersect()` means it is not necessary to use `order` afterwards and makes things easier. Hadn't thought of that... Thanks! – Robert Hacken Sep 15 '22 at 17:42
I found the suggestions of Robert Hacken very helpful, but unfortunately it goes wrong again when not all locations are found in the same polygon. For example when adding Hollum (found on a island in the North) with lat=53.440278 and long=5.639722. – WJH Sep 20 '22 at 07:54
Is there any solution for this? – WJH Sep 20 '22 at 08:00
1

Yes, there is, see the end of my updated answer. The problem was that the island point was outside the area of Netherlands returned by `ne_countries()`, so matching by `st_intersection()` was not succesful. – Robert Hacken Sep 20 '22 at 17:58

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

1 Answers1