6.1 Data set: wrangling

library(sp)
library(SpatialEpi)
library(sf)
library(leaflet)

The data set cover Scoltland counties:

• Observed and expected cases of lip cancer in males

• % of opop in Agriculture, Fishing, Forestry (AFF)

6.1.1 Map of Scotland counties

We will be using both {sp} and {sf} just for the sake of understanding:

data(scotland)
map <- scotland$spatial.polygon
plot(map)

The workflow in {sp} an start with getting all EPSG codes and proj4 fron rgdal:

# sp
codes <- rgdal::make_EPSG()
codes[which(codes$code == "27700"), ]

##       code                           note
## 6267 27700 OSGB36 / British National Grid
##                                                                                                                  prj4
## 6267 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
##               prj_method
## 6267 Transverse Mercator

codes[codes$code == "4326",]

##      code   note                                          prj4 prj_method
## 2180 4326 WGS 84 +proj=longlat +datum=WGS84 +no_defs +type=crs     (null)

# saving it for other purposes
map_3 <- map_2 <-map
# change in unit was done by "hand"
proj4string(map) <- "+proj=tmerc +lat_0=49 +lon_0=-2
+k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36
+units=km +no_defs"

## Warning in ReplProj4string(obj, CRS(value)): A new CRS was assigned to an object with an existing CRS:
## +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransform

An other option if you want to stay in {sp} but avoid {rgdal}. You have to keep in mind that {sp} use S4 classes.

# partially sf
map_2 <- map
(proj4string(map_2) <- sp::CRS(SRS_string = sf::st_crs(27700)$wkt))

## Warning in `proj4string<-`(`*tmp*`, value = new("CRS", projargs = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs")): A new CRS was assigned to an object with an existing CRS:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=km +no_defs
## without reprojecting.
## For reprojection, use function spTransform

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["OSGB36 / British National Grid",
##     BASEGEOGCRS["OSGB36",
##         DATUM["Ordnance Survey of Great Britain 1936",
##             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4277]],
##     CONVERSION["British National Grid",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",49,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-2,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996012717,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",400000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",-100000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
##         BBOX[49.75,-9,61.01,2.01]],
##     ID["EPSG",27700]]

Doing it with full {sf}:

# More sf 
map_3 <- st_as_sf(map)
map_3 <- st_set_crs(map_3, 27700) # print st_set_crs()

## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that

map_3$name <- sapply(map@polygons, function(x) x@ID)
# an other way of doing it provided by Paula Moraga
# sapply(slot(map, "polygons"), function(x){slot(x, "ID")})
# I did not know about slot()

More on proj4string and crs as wkt here: https://r-spatial.org/r/2020/03/17/wkt.html

For leaflet we will need to transform the CRS:

# require rgdal
map <- sp::spTransform(map,
                   CRS("+proj=longlat +datum=WGS84 +no_defs"))

## Warning: PROJ support is provided by the sf and terra packages among others

map_3 <- sf::st_transform(map_3, 4326)

6.1.2 Data

With {sp}:

d <- scotland$data[,c("county.names", "cases", "expected", "AFF")]
names(d) <- c("county", "Y", "E", "AFF")
d$SIR <- d$Y / d$E
rownames(d) <- d$county
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE)

With {sf}:

# you can also use dplyr verbs
map_3m <- merge(map_3, d, by.x = "name", by.y = "county")