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")
How do I transfer the meuse.sf['zinc']
values to the Voronoi and visualize the polygons colored like the points please?
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