Posted to Tutorials  /  Tags: , ,

Voronoi Diagram and Delaunay Triangulation in R

The deldir package by Rolf Turner makes the calculations and plotting straightforward, with a few lines of code.

delaunay-featured

Read more on Wikipedia.A Voronoi diagram splits up a plane based on a set of original points. Each polygon, or Voronoi cell, contains an original point and all areas that are closer to that point than any other. The Delaunay triangulation of the same set of original points is the dual graph of the Voronoi diagram.

We see them around these parts occasionally.

The main thing you need to know for this R tutorial though is that the deldir package by Rolf Turner does most of the work for you. Calculations and plotting only takes a couple lines of code.

Setup

R is an open source statistical computing language useful for analysis. You can also use it as a visualization tool.This tutorial makes use of R and a handful of packages, so install these first if you haven’t already.

To install R, go to the homepage, click on the link to download R, and select the file that corresponds to your system. It’s typically a straightforward installation with only a few clickety clicks.

With R installed, open the console (i.e. the app) and install the necessary packages using the install.packages() function. You definitely need the deldir package:

install.packages("deldir", dependencies=TRUE)

The second part of this tutorial does some stuff with spatial data and you’ll need some packages for that.

install.packages(c("sp", "maps", "maptools", "mapproj"), dependencies=TRUE)

If you don’t care about applying the Voronoi and Delaunay to spatial data though, you can skip these packages if you like.

Unzip the file and set your working directory in R to the resulting directory.Finally, download the source for this tutorial, and use it to follow along. The data also includes data that you will plot later on.

Calculations and Plotting

Start with some random data on an x-y plane.

# Generate points
x <- rnorm(500, 0, 1.5)
y <- rnorm(500, 0, 1)

This gives you 500 x-coordinates and 500 y-coordinates sampled from normal distributions with mean 0 and standard deviation 1.5, and mean 0 and standard deviation 1, respectively.

To calculate the Voronoi tessellation and the Delaunay triangulation, just call the deldir() function on x and y.

# Calculate tessellation and triangulation
vtess <- deldir(x, y)

Then plot on x and y to setup a plotting window, without any points.

plot(x, y, type="n", asp=1)

You should see a blank plot space. Add the points:

points(x, y, pch=20, col="red", cex=0.5)

Then plot the results of deldir(), specifying wlines as “tess” to to plot a Voronoi diagram and add to TRUE too add the diagram to the current plotting window.

plot(vtess, wlines="tess", wpoints="none", number=FALSE, add=TRUE, lty=1)

The result:

Voronoi diagram

Nice.

Similarly, you can plot the triangulation without any new calculations. Just set wlines to “triang” instead of “tess” like in the previous function call.

# Delaunay triangulation
plot(x, y, type="n", asp=1)
plot(vtess, wlines="triang", wpoints="none", number=FALSE, add=TRUE, lty=1)
points(x, y, pch=20, col="black", cex=0.5)

The result:

Delaunay triangulation

The process isn’t much different with real data. Let’s try it out with airport locations in the conterminous United States like in Mike Bostock’s d3.js example.

Be sure to set your working directory to where you saved the tutorial download.First, load the data with read.csv().

# Load data
airports <- read.csv("data/airport-locations.tsv", sep="\t", stringsAsFactors=FALSE)

Sorry non-conterminous United States.This dataset includes all US locations, which means it includes Hawaii, Alaska, and Puerto Rico. For this example, we’re only interested in the conterminous United States, so filter accordingly.

source("lib/latlong2state.R")
airports$state <- latlong2state(airports[,c("longitude", "latitude")])
airports_contig <- na.omit(airports)

I covered this in more detail in the gridded maps tutorial.The latlong2state() function, which came directly from a StackOverflow answer. it takes a data frame with longitude and latitude, in that order, and returns the the state that a longitude-latitude point is in. If the location is not in the conterminous United States the function returns NA. So run latlong2state() on the airports data frame, and then omit any rows that have NA state.

This gives you a data frame that looks like this:

   latitude  longitude       state
1  31.95376  -89.23450 mississippi
2  30.68586  -95.01793       texas
3  38.94575 -104.56989    colorado
4  42.74135  -78.05208    new york
5  30.68801  -81.90594     florida
6  34.49167  -88.20111 mississippi
7  32.85049  -86.61145     alabama
8  43.08751  -88.17787   wisconsin
...

We’ve looked at map projections quite a few times in previous tutorials. You can find more here in the geographic paths tutorial. A projection isn’t necessary for this example, but we’re trying to duplicate the Bostock example. Plus, it looks nice.There are 2,988 locations total. Project the points to Albers:

# Projection
library(mapproj)
airports_projected <- mapproject(airports_contig$longitude, airports_contig$latitude, "albers", param=c(39,45))

Like before, you can plot the (projected) points.

# Plot
par(mar=c(0,0,0,0))
plot(airports_projected, asp=1, type="n", bty="n", xlab="", ylab="", axes=FALSE)
points(airports_projected, pch=20, cex=0.1, col="red")

Think of the data as a standard x-y data frame. It just happens to be geographic, where x is longitude and y is latitude.

Point map

With this in mind, you calculate the Voronoi tessellation and Delaunay triangulation in the same way with deldir().

# Voronoi
vtess <- deldir(airports_projected$x, airports_projected$y)
plot(vtess, wlines="tess", wpoints="none", number=FALSE, add=TRUE, lty=1)

Note that add in deldir() is set to TRUE to add to the current plot (a map in this case). You get what you expect.

Voronoi airports

Same deal with the Delaunay triangulation.

# Delaunay triangulation
vtess <- deldir(airports_projected$x, airports_projected$y)
plot(vtess, wlines="triang", wpoints="none", number=FALSE, add=TRUE, lty=1)

Yay.

Delaunay triangulation airports

Wrapping up

If you need a Voronoi tessellation or a Delaunay triangulation based on a set of points in a two-dimensional space, the deldir package has you covered in a couple of lines of code. The package provides other functions too like plot.tile.list() which lets you color cells on a per-cell basis. Get all the details on the package page.

You might also be interested in Paul Murrell’s implementation of Voronoi treemaps.

About the Author

Nathan Yau is a statistician who works primarily with visualization. He earned his PhD in statistics from UCLA, is the author of two best-selling books — Data Points and Visualize This — and runs FlowingData. Introvert. Likes food. Likes beer. Follow him @flowingdata.

Become a member. Learn to visualize your data. Support FlowingData.

Join Now

Membership

This is for people who want to learn to make and design data graphics. Your support goes directly to FlowingData, an independently run site.

What You Get

  • Instant access to tutorials on how to make and design data graphics
  • Source code and files to use with your own data
  • Four-week course on visualization in R
  • Hand-picked links and resources from around the web

More Tutorials See All →

Small Multiples in R

Make a lot of charts at once, line them up in a grid, and you can make quick comparisons across several categories.

How to Download and Use Online Data with Arduino

Before you can do anything with data, you have to get it into the application. Working with an Arduino is no different. Although the process is changes, if you’re used to working with desktop software.

R Cheat Sheet and Guide for Graphical Parameters

You can customize graphics in R with par(), but the docs are mostly text and just organized alphabetically. Here is a more visual reference, categorized by what you can change. Plus, a one-page printout.

Choropleth Maps and Shapefiles in R

Fill those empty polygons with color, based on shapefile or external data.