John Snow’s Cholera Map

References

Prepare dataset from Robin’s Blog

## Download the zip archive of geographical information
download.file(url      = "http://www.rtwilson.com/downloads/SnowGIS_v2.zip",
              destfile = "SnowGIS_v2.zip")

## Unzip
unzip(zipfile = "SnowGIS_v2.zip")

## List files in the unzipped folder
dir(path = "./SnowGIS")
 [1] "Cholera_Deaths.dbf"          "Cholera_Deaths.prj"          "Cholera_Deaths.sbn"         
 [4] "Cholera_Deaths.sbx"          "Cholera_Deaths.shp"          "Cholera_Deaths.shx"         
 [7] "OSMap_Grayscale.tfw"         "OSMap_Grayscale.tif"         "OSMap_Grayscale.tif.aux.xml"
[10] "OSMap_Grayscale.tif.ovr"     "OSMap.tfw"                   "OSMap.tif"                  
[13] "Pumps.dbf"                   "Pumps.prj"                   "Pumps.sbx"                  
[16] "Pumps.shp"                   "Pumps.shx"                   "README.txt"                 
[19] "SnowMap.tfw"                 "SnowMap.tif"                 "SnowMap.tif.aux.xml"        
[22] "SnowMap.tif.ovr"            

## Load rgdal package
library(rgdal)

## Check tiff file
GDALinfo("./SnowGIS/OSMap.tif")
rows        1053 
columns     1169 
bands       1 
lower left origin.x        528765 
lower left origin.y        180466 
res.x       1 
res.y       1 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m
+no_defs 
file        ./SnowGIS/OSMap.tif 
apparent band summary:
  GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1   Byte           TRUE         255        128        128
apparent band statistics:
  Bmin Bmax Bmean Bsd
1    0  255    NA  NA
Metadata:
AREA_OR_POINT=Area 

## Load gray scale image
OSMap <- readGDAL("./SnowGIS/OSMap_Grayscale.tif")
./SnowGIS/OSMap_Grayscale.tif has GDAL driver GTiff 
and has 1070 rows and 1169 columns

## Set margins to zero
par(mar = c(0,0,0,0))

## Plot it
image(OSMap, col = grey(1:999/1000))

## Load maptools package
library(maptools)

## Get shapefile header information (cholera death sites)
getinfo.shape("./SnowGIS/Cholera_Deaths.shp")
Shapefile type: Point, (1), # of Shapes: 250

## Load geographical information of sites of cholera deaths
Cholera_Deaths <- readShapePoints(fn = "./SnowGIS/Cholera_Deaths.shp")

## Plot sites of cholera deaths
plot(Cholera_Deaths, pch = 16, cex = log(Cholera_Deaths$Count), col = "red", add = TRUE)

## Get shapefile header information (pump sites)
getinfo.shape("./SnowGIS/Pumps.shp")
Shapefile type: Point, (1), # of Shapes: 8

## Load geographical information of sites of pumps
Pumps <- readShapePoints(fn = "./SnowGIS/Pumps.shp")

## Plot sites of pumps
plot(Pumps, add = TRUE, pch = 16, col = "blue")

plot of chunk unnamed-chunk-3