Getting data into R

Data import from T:/ drive

John Snow GIS data saved on T:/ drive along with a nice ‘mapping cholera’ exercise. For now I’m only grabbing the GIS components. This was all prepared by

library(sf)
library(raster)

tt <- 'T:/epiprojs/KRAMER_GIS/JohnSnowCholera/SnowGIS/'

dir(tt)
##  [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.tif.aux.xml"           "OSMap.tif.ovr"              
## [11] "OSMap_Grayscale.tfw"         "OSMap_Grayscale.tif"        
## [13] "OSMap_Grayscale.tif.aux.xml" "OSMap_Grayscale.tif.ovr"    
## [15] "Pumps.dbf"                   "Pumps.prj"                  
## [17] "Pumps.sbx"                   "Pumps.shp"                  
## [19] "Pumps.shx"                   "README.txt"                 
## [21] "SnowMap.tfw"                 "SnowMap.tif"                
## [23] "SnowMap.tif.aux.xml"         "SnowMap.tif.ovr"

So there are 2 shapefiles:

  • Cholera_Deaths.shp
  • Pumps.shp

And there are several raster files:

  • OSMap.tif
  • OSMap_Grayscale.tif
  • SnowMap.tif

Read in files

**Vector data:**

deaths <- st_read(paste(tt, 'Cholera_Deaths.shp', sep = ''))
## Reading layer `Cholera_Deaths' from data source `T:\epiprojs\KRAMER_GIS\JohnSnowCholera\SnowGIS\Cholera_Deaths.shp' using driver `ESRI Shapefile'
## Simple feature collection with 250 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 529160.3 ymin: 180857.9 xmax: 529655.9 ymax: 181306.2
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
pumps <- st_read(paste(tt, 'Pumps.shp', sep =''))
## Reading layer `Pumps' from data source `T:\epiprojs\KRAMER_GIS\JohnSnowCholera\SnowGIS\Pumps.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 529183.7 ymin: 180660.5 xmax: 529748.9 ymax: 181193.7
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs

**Raster data:**

NOTE: the SnowMap.tif is 26 MB whereas the other ‘recent’ maps are <.5 MB

OSMap <- raster(paste(tt, 'OSMap.tif', sep =''))
OSMap_gs <- raster(paste(tt, 'OSMap_Grayscale.tif', sep = ''))
snow <- raster(paste(tt, 'SnowMap.tif', sep = ''))

Plot maps of data

First, plotting over modern day London map (in color and greyscale)

library(tmap)

p <- tm_shape(pumps) + tm_dots(col = 'blue', size = 1, shape = 22, legend.show = F)

d <- tm_shape(deaths) + tm_bubbles(size = 'Count', col = 'red')

bm <- tm_shape(OSMap) + tm_raster()

bm_gs <- tm_shape(OSMap_gs) + tm_raster()


bm + d + p

bm_gs + d + p

Now plotting Snow’s original alone, and then with data points

tm_shape(snow) + tm_raster(palette =  '-Greys')

s <- tm_shape(snow) + tm_raster(palette = '-Greys')

s + d + p

Saving data to Course folder

st_write(deaths, 'CholeraDeaths.gpkg')
st_write(pumps, 'Pumps.gpkg')
writeRaster(OSMap, 'OSMap_color.tif')
writeRaster(OSMap_gs, 'OSMap_gs.tif')
writeRaster(snow, 'JS_map.tif')

Other Cholera data

NOTE: The below didn’t end up working because it is not clear what projection the cholera data are in, therefore I can’t make it align with the other data!

In addition to the above GIS data, there is an R package dedicated to analysis of the Snow cholera data: https://github.com/lindbrook/cholera

From this package there are veroni neighborhoods which I’m grabbing to add a polygon to the collection for lab exercise…

library(cholera)

x <- neighborhoodVoronoi() # this produces a bunch of stuff including 'coordinates'
x2 <- x$coordinates # this is a list of polygon vertice coordiantes

Here I created a voroni neighborhood, the vertices of which are in the coordinates object within the x object. To make this into an sf polygon I do need lists of points, but they need to be matrix not data.frame as they are here. So converting them…

x3 <- lapply(x2, function(x) {if(any(class(x)=="data.frame")) as.matrix(x) else x})

Now creating a polygon file: (NOTE DISABLED FOR NOW)

x3_p <- st_polygon(x3)
x3_p <- st_sf(geom = st_sfc(x3_p))
class(x3_p)