John Snow, and OpenStreetMap

While I was working for a training on data visualization, I wanted to get a nice visual for John Snow’s cholera dataset. This dataset can actually be found in a great package of famous historical datasets.


One can easily visualize the deaths, on a simplified map, with the streets (here simple grey segments, see Vincent Arel-Bundock’s post)

plot(Snow.deaths[,c("x","y")], col="red", pch=19, cex=.7,xlab="", ylab="", xlim=c(3,20), ylim=c(3,20))
slist <- split(Snow.streets[,c("x","y")],as.factor(Snow.streets[,"street"]))
invisible(lapply(slist, lines, col="grey"))

Of course, one might add isodensity curves (estimated using kernels)

kde2d <- bkde2D(Snow.deaths[,2:3], bandwidth=c(0.5,0.5))
contour(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat, add=TRUE)

Now, what if we want to visualize that dataset on a nice background, from Google Maps, or OpenStreetMaps? The problem here is that locations are in a weird coordinate representation system. So let us use a different dataset. For instance, on Robin Wilson’s blog, one can get datasets in a more traditional representation (here the epsg 27700). We can extract the dataset from

> library(maptools)
> setwd("/cholera/")
> deaths <- readShapePoints("Cholera_Deaths")
> head(deaths@coords)
coords.x1 coords.x2
0  529308.7  181031.4
1  529312.2  181025.2
2  529314.4  181020.3
3  529317.4  181014.3
4  529320.7  181007.9
5  529336.7  181006.0

Then, we need our background,

map = openmap(c(lat= 51.516,   lon= -.141),
              c(lat= 51.511,   lon= -.133))
map=openproj(map, projection = "+init=epsg:27700") 
points(deaths@coords,col="red", pch=19, cex=.7 )

If we zoom in (the code above will be just fine), we get

And then, we can compute the density

kde2d <- bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])))

based on the same function as before (here I use marginal cross-validation techniques to get optimal bandwidths). To get a nice gradient, we can use

clrs=colorRampPalette(c(rgb(0,0,1,0), rgb(0,0,1,1)), alpha = TRUE)(20)

and finally, we add it on the map

image(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat, add=TRUE,col=clrs)
contour(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat, add=TRUE)

Cite this blog post
Arthur Charpentier (2015, February 27). John Snow, and OpenStreetMap. Freakonometrics. Retrieved May 23, 2024, from

11 thoughts on “John Snow, and OpenStreetMap”

  1. Hi,
    I try to practice these steps, everything is perfect but there is an error:
    > points(deaths@coords, col=”red”, pch=19, cex=0.7)
    Error in points(deaths@coords, col = “red”, pch = 19, cex = 0.7) :
    trying to get slot “coords” from an object (class “data.frame”) that is not an S4 object

    I doubt that this error is related to the new version of dataset :, could you please share the dataset which you use or guide me the right direction?


Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.