## 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")