In the HistData package for R (https://r-forge.r-project.org/R/?group_id=574) I have the data sets related to John Snow's map of the cholera outbreak in London, 1854. I believe they are authoritative, having been carefully digitized under supervision of Walter Tobler. Some details on these data sets are described by John Mackenzie, at http://www1.udel.edu/johnmack/frec480/cholera/cholera2.html.

Unfortunately, the coordinates of deaths, pumps, and streets use an arbitrary coordinate system, not map coordinates suitable for other GIS applications or mapping software in R (spatial packages, ggmap, etc.)

In http://freakonometrics.hypotheses.org/19213 Arthur Charpentier uses ggmap with a version of the John Snow data from http://www.rtwilson.com/downloads/SnowGIS_v2.zip. The Cholera_Deaths.shp file, however lists only 489 deaths, not the 578 I have recorded in HistData::Snow.deaths.

One idea is to find the relations between the means and standard deviations of the (x,y) coordinates and rescale linearly, but maybe there is a better way?

Here's what I've tried so far

> data(Snow.deaths, package="HistData")
> D <- Snow.deaths[,2:3]
> colMeans(D)
       x        y 
13.03312 11.69721 
> var(D)
          x         y
x 3.8150987 0.3802654
y 0.3802654 2.7213828

Read the Cholera_deaths file

> folder <- "C:/Dropbox/R/data/Snow/SnowGIS_v2/SnowGIS"
> library(maptools)
> deaths <- readShapePoints(file.path(folder, "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
> # deaths has only 250 observations; 489 deaths
> sum(deaths@data$Count)
[1] 489

 > # try to relate to Snow.deaths
> X <- deaths@coords
> colnames(X) <- c("x", "y")
> XX <- data.frame(X, Freq=deaths@data$Count)
> XX <- vcdExtra::expand.dft(XX)
> colMeans(XX)
       x        y 
529414.8 181031.9 
> var(XX)
          x        y
x 10813.816 1521.693
y  1521.693 6227.924

OK, then I try to re-scale D to have the same means and standard deviations as XX, but something doesn't work correctly here-- the column means of Dscaled should have turned out to be equal to those of XX:

> # scale D to have the same means and standard deviations as XX
> Dscaled <- scale(D, center=TRUE, scale=TRUE)
> Dscaled <- scale(Dscaled, center=colMeans(XX), scale=sqrt(diag(var(XX))))
> colMeans(Dscaled)
        x         y 
-5091.040 -2293.947 

EDIT: It might be helpful in this problem to see Snow's map as drawn by the new function, SnowMap(axis.labels=TRUE) now in the development version of HistData (rev 102) on R-Forge. The axis labels show the origin of the coordinate system in the lower left corner as they are in my data Snow.* data sets.


  • I've just been trying to scale the pumps from each data set to match up. I don't believe the line in help(Snow.pumps) about the coordinates being 100 metres, since a scale of about 54 (and a translation) do the best job of mapping these to the shapefile's UK grid coordinates (which definitely are in metres). Even then the points don't overlap exactly, some other rotation/skew is clearly present. Since there are fewer pumps its possible to identify corresponding pumps in each data set and compute the shift/translate for those. – Spacedman Dec 24 '16 at 8:33
  • I assume you've looked at HistData/inst/doc/Snow_deaths-duplicates.html and found it unhelpful? – user1462 Dec 24 '16 at 12:55
  • It also occurred to me that I could obtain the linear transformation of coordinates in my Snow.* files to those in a GIS-based map with the locations of two pumps, or three to check the accuracy. Unfortunately, there are no labels for the pumps in the SnowGIS files, and I haven't seen an example of how to plot them so that I can compare them visually. – user101089 Dec 24 '16 at 19:56
  • 1
    For a second there after reading your title, i thought you wanted to map coordinates in Westeros. – user35594 Dec 24 '16 at 20:16

Perhaps evaluate the shapefile from http://donboyes.com/2011/10/14/john-snow-and-serendipity it has 578 points.

I don't think trying to relate the HistData Snow Deaths to Robin Wilson's (@robintw) version will work as the shapefile contains a single point coordinate for multiple deaths at a single address, instead of the multiple points stacked back from the street in the map.

Robin's version is definitely missing a lot of points. From a quick look, there's quite a few outlying single deaths missed. Another problem is closer to the centre of the map where it was not properly edge matched when put together (this is also visible in the Wikipedia map) and this obscures a number of points.

Extract of map supplied in the download:

enter image description here

Extract from UCLA version:

enter image description here

| improve this answer | |

To complete the answer to this question, the following code finds the linear transformation of the coordinates in the original Tobler files (in HistData) and those proived by Don Boyes.

folder <- "C:/Dropbox/R/data/Snow/snow_shp"
deaths <- readShapePoints(file.path(folder, "deaths_gcs"))
data(Snow.deaths, package="HistData")
X <- deaths@coords
D <- Snow.deaths[,2:3]

Then, correlate and regress D[,1] on X[,1] and D[,2] on X[,2]. The linear transformation is given by the regression coefficients.

> cor(D[,1], X[,1])
[1] 0.9999664
> cor(D[,2], X[,2])
[1] 0.9995559
> # linear transformations to GIS coords
> lm(D[,1] ~ X[,1])

lm(formula = D[, 1] ~ X[, 1])

(Intercept)       X[, 1]  
      185.4       1264.7  

> lm(D[,2] ~ X[,2])

lm(formula = D[, 2] ~ X[, 2])

(Intercept)       X[, 2]  
    -105441         2047  
| improve this answer | |

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.