7

For a point data set with temperature data I have, I want to create their Thiessen Polygons in R. So far, the only algorithm I have found is that of Carson Farmer (http://carsonfarmer.com/2009/09/voronoi-polygons-with-r/). I would like to know if there is a ready to use function that creates Thiessn Polygons?

5

Thiessen polygons are Voronoi diagrams - there is a 'voronoi' package available in the CRAN archives (not the main repository), but the 'deldir' package does the same job.

require(deldir)

# Create some points
x <- c(32.5, 32.1, 33.5, 32.2, 33.0)
y <- c(-2.2, -3.3, -2.3, -2.9, -3.0)

# Calculate the Delaunay triangulation, then the tiles.
z <- deldir(x,y,rw=c(32.0,33.6,-3.4,-2.1))
w <- tile.list(z)

# Make a list of pretty colours, and use 'em to plot:
wcols <- topo.colors(5) 
plot(w, fillcol=wcols, close=TRUE)

enter image description here

Now, the 'w' object is a 'tile list', not an sp object, but the tiles could be turned into a spatial object pretty easily, using the x / y components of the list:

> str(w[1])
List of 1
 $ :List of 5
  ..$ ptNum: int 1
  ..$ pt   : Named num [1:2] 32.5 -2.2
  .. ..- attr(*, "names")= chr [1:2] "x" "y"
  ..$ x    : num [1:5] 33 32 32 32.6 33
  ..$ y    : num [1:5] -2.1 -2.1 -2.4 -2.67 -2.46
  ..$ bp   : logi [1:5] TRUE TRUE TRUE FALSE FALSE
  • I had a script saved to create an sp object using a tesselation from the spatstat package. See here for your example points. I'm pretty sure spatstat just uses the same deldir library, but I'm lazy and couldn't figure out how to edit the functions to get them to work directly with the tile list from deldir in a quick attempt. – Andy W Feb 24 '15 at 13:17
  • Thanks for the script - if I can get around to making it work with deldir's output I'll update here. – Simbamangu Mar 3 '15 at 14:02
4

Alternatively, you could use dismo package. It accepts SpatialPoints* or matrix as input data.

library('dismo')
x <- c(32.5, 32.1, 33.5, 32.2, 33.0)
y <- c(-2.2, -3.3, -2.3, -2.9, -3.0)
points <- matrix(c(x,y), ncol=2)
vor <- voronoi(points)
spplot(vor, "id")

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.