rplotcolorspolygonvoronoi

Assign color to Voronoi


I would like to assign each Thiessen polygon the same color as the point used to generate the tessellation. We can use the meuse dataset.

library(sf) # 'simple features' representations of spatial objects
#- read
data("meuse", package = "sp")
#- to sf
meuse.sf <- st_as_sf(meuse, coords = c("x","y"))
#- crs
st_crs(meuse.sf) <- 28992

#- Voronoi tesselation
pnts <- st_union(meuse.sf)
voronoi_grid <- st_voronoi(pnts)

#-- plot
plot(voronoi_grid, col = NA, #meuse.sf[,"zinc"],
      xlim=c(178605,181390),
      ylim=c(329714,333611))
points(st_coordinates(meuse.sf)[,1], st_coordinates(meuse.sf)[,2],
       pch=21,
       bg=sp::bpy.colors(length(meuse.sf$zinc))[rank(meuse.sf$zinc)],
       cex=0.4*meuse.sf$zinc/max(meuse.sf$zinc))
#grid()
title("Tesselated Voronoi Surface")

enter image description here

How do I transfer the meuse.sf['zinc'] values to the Voronoi and visualize the polygons colored like the points please?


Solution

  • You can extract polygons from the GEOMETRYCOLLECTION that's returned by st_voronoi(), convert the resulting set to an sf object and spatially join it to (a subset of) meuse.sf:

    library(sf) 
    data("meuse", package = "sp")
    meuse.sf <- st_as_sf(meuse, coords = c("x","y"))
    st_crs(meuse.sf) <- 28992
    
    pnts <- st_union(meuse.sf)
    voronoi_grid <- st_voronoi(pnts)
    
    # check returned object
    voronoi_grid
    #> Geometry set for 1 feature 
    #> Geometry type: GEOMETRYCOLLECTION
    #> Dimension:     XY
    #> Bounding box:  xmin: 174708 ymin: 325817 xmax: 185287 ymax: 337508
    #> Projected CRS: Amersfoort / RD New
    #> GEOMETRYCOLLECTION (POLYGON ((174708 332989.6, ...
    
    # extract POLYGONs, turn resulting geomerty set to sf and 
    # spatially join to subset of meuse.sf:
    voronoi_sf <- st_sf(geometry = st_collection_extract(voronoi_grid, "POLYGON")) |>
      st_join(meuse.sf[,"zinc"])
    voronoi_sf
    #> Simple feature collection with 155 features and 1 field
    #> Geometry type: POLYGON
    #> Dimension:     XY
    #> Bounding box:  xmin: 174708 ymin: 325817 xmax: 185287 ymax: 337508
    #> Projected CRS: Amersfoort / RD New
    #> First 10 features:
    #>    zinc                       geometry
    #> 1   553 POLYGON ((174708 332989.6, ...
    #> 2   332 POLYGON ((178621.5 330091.1...
    #> 3   907 POLYGON ((174708 335077.5, ...
    #> 4   577 POLYGON ((174708 333882.1, ...
    #> 5   560 POLYGON ((178796.1 330570.9...
    #> 6   601 POLYGON ((177424 325817, 17...
    #> 7   783 POLYGON ((174708 325817, 17...
    #> 8   342 POLYGON ((178775.3 330076.5...
    #> 9  1383 POLYGON ((178463.4 331081.4...
    #> 10  420 POLYGON ((178770.5 330462.1...
    
    plot(voronoi_sf[,"zinc"],
         xlim=c(178605,181390),
         ylim=c(329714,333611))
    

    Created on 2023-09-05 with reprex v2.0.2