For this assignment you will use the library(HistData) and the dataset Snow to recreate the famed John Snow map using the ggmap() package. Feel free to embelish and annotate the map to any degree you wish. Once completed upload the Rmarkdown document to Rpubs.com (accounts are free) and obtain a “share” link. The address of that link should be submitted on Moodle.
## 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"
## [3] "Cholera_Deaths.sbn" "Cholera_Deaths.sbx"
## [5] "Cholera_Deaths.shp" "Cholera_Deaths.shx"
## [7] "OSMap.tfw" "OSMap.tif"
## [9] "OSMap_Grayscale.tfw" "OSMap_Grayscale.tif"
## [11] "OSMap_Grayscale.tif.aux.xml" "OSMap_Grayscale.tif.ovr"
## [13] "Pumps.dbf" "Pumps.prj"
## [15] "Pumps.sbx" "Pumps.shp"
## [17] "Pumps.shx" "README.txt"
## [19] "SnowMap.tfw" "SnowMap.tif"
## [21] "SnowMap.tif.aux.xml" "SnowMap.tif.ovr"
library(ggmap)
library(maptools)
library(rgdal)
library(maptools)
library(leaflet)
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
OSMap <- readGDAL("./SnowGIS/OSMap_Grayscale.tif")
## ./SnowGIS/OSMap_Grayscale.tif has GDAL driver GTiff
## and has 1070 rows and 1169 columns
par(mar = c(0,0,0,0))
image(OSMap, col = grey(1:999/1000))
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")
london_map <- get_map(c(-.137,51.513), zoom=16)
london <- ggmap(london_map)
setwd("C:/Users/Vinay/SnowGIS")
deaths <- readShapePoints("Cholera_Deaths")
pumps <- readShapePoints("Pumps")
df_deaths <- data.frame(deaths@coords)
df_pumps <- data.frame(pumps@coords)
tmp <- rbind(df_deaths, df_pumps)
tmp$type <- c(rep('death', times=dim(df_deaths)[1]),
rep('pump', times=dim(df_pumps)[1]))
coordinates(df_deaths)=~coords.x1+coords.x2
proj4string(df_deaths)=CRS("+init=epsg:27700")
df_deaths = spTransform(df_deaths,CRS("+proj=longlat +datum=WGS84"))
df=data.frame(df_deaths@coords)
lng=df$coords.x1
lat=df$coords.x2
coordinates(tmp)=~coords.x1+coords.x2
proj4string(tmp)=CRS("+init=epsg:27700")
tmp = spTransform(tmp, CRS("+proj=longlat +datum=WGS84"))
tmp <- data.frame(tmp@coords, type=tmp@data$type)
snow.plot <- london
snow.plot
london +
geom_point(mapping=
aes(x=coords.x1, y=coords.x2, col=type),
data=tmp)
snow.plot2 <- snow.plot + geom_density2d(data = tmp[tmp$type ==
"death", ], aes(x = coords.x1, y = coords.x2),
size = 0.3)
snow.plot2
snow.plot3 <- snow.plot2 + stat_density2d(data = tmp[tmp$type ==
"death", ], aes(x = coords.x1, y = coords.x2,
fill = ..level.., alpha = ..level..),
size = 0.01, bins = 16, geom = "polygon") +
scale_fill_gradient(low = "red",
high = "red", guide = FALSE) +
scale_alpha(range = c(0, 0.3), guide = FALSE)
snow.plot3