library(tidyverse)
library(mdsr)
library(sf)
plot(CholeraDeaths["Count"])
Spatial Data
The new version of the book is now available mdsr2e. The new chapter Chapter 17 Working with geospatial data used the new R package sf and could be further updated to the the new R package cholera.
It is probably best to do this in an R Project.
download.file("http://rtwilson.com/downloads/SnowGIS_SHP.zip",
dest="SnowGIS.zip", mode="wb")
unzip("SnowGIS.zip")
dsn <- fs::path("SnowGIS_SHP")
list.files(dsn) [1] "Cholera_Deaths.dbf" "Cholera_Deaths.prj"
[3] "Cholera_Deaths.sbn" "Cholera_Deaths.sbx"
[5] "Cholera_Deaths.shp" "Cholera_Deaths.shx"
[7] "OSMap_Grayscale.tfw" "OSMap_Grayscale.tif"
[9] "OSMap_Grayscale.tif.aux.xml" "OSMap_Grayscale.tif.ovr"
[11] "OSMap.tfw" "OSMap.tif"
[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" st_layers(dsn)Driver: ESRI Shapefile
Available layers:
layer_name geometry_type features fields crs_name
1 Pumps Point 8 1 OSGB36 / British National Grid
2 Cholera_Deaths Point 250 2 OSGB36 / British National GridCholeraDeaths <- st_read(dsn, layer = "Cholera_Deaths")Reading layer `Cholera_Deaths' from data source
`/home/esuess/classes/2023-2024/Fall 2023/Stat651/Presentations/08_spatial/spatial4/SnowGIS_SHP'
using driver `ESRI Shapefile'
Simple feature collection with 250 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 529160.3 ymin: 180857.9 xmax: 529655.9 ymax: 181306.2
Projected CRS: OSGB36 / British National Gridclass(CholeraDeaths)[1] "sf" "data.frame"CholeraDeathsSimple feature collection with 250 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 529160.3 ymin: 180857.9 xmax: 529655.9 ymax: 181306.2
Projected CRS: OSGB36 / British National Grid
First 10 features:
Id Count geometry
1 0 3 POINT (529308.7 181031.4)
2 0 2 POINT (529312.2 181025.2)
3 0 1 POINT (529314.4 181020.3)
4 0 1 POINT (529317.4 181014.3)
5 0 4 POINT (529320.7 181007.9)
6 0 2 POINT (529336.7 181006)
7 0 2 POINT (529290.1 181024.4)
8 0 2 POINT (529301 181021.2)
9 0 3 POINT (529285 181020.2)
10 0 2 POINT (529288.4 181031.8)summary(CholeraDeaths) Id Count geometry
Min. :0 Min. : 1.000 POINT :250
1st Qu.:0 1st Qu.: 1.000 epsg:27700 : 0
Median :0 Median : 1.000 +proj=tmer...: 0
Mean :0 Mean : 1.956
3rd Qu.:0 3rd Qu.: 2.000
Max. :0 Max. :15.000 ggplot(CholeraDeaths) +
geom_sf()
In the new edition of the book the authors use Open Street Maps instead of Google Maps API.
library(ggspatial)
library(prettymapr)
ggplot(CholeraDeaths) +
annotation_map_tile(type = "osm", zoomin = 0) +
geom_sf(aes(size = Count), alpha = 0.7)Loading required namespace: rasterZoom: 17
st_bbox(CholeraDeaths) xmin ymin xmax ymax
529160.3 180857.9 529655.9 181306.2 Projections
library(mapproj)Loading required package: maps
Attaching package: 'maps'The following object is masked from 'package:purrr':
maplibrary(maps)
map("world", projection = "mercator", wrap = TRUE)Warning in map("world", projection = "mercator", wrap = TRUE): projection
failed for some data
map("world", projection = "cylequalarea", param = 45, wrap = TRUE)
map(
"state", projection = "lambert",
parameters = c(lat0 = 20, lat1 = 50), wrap = TRUE
)
map(
"state", projection = "albers",
parameters = c(lat0 = 20, lat1 = 50), wrap = TRUE
)
st_crs(CholeraDeaths)Coordinate Reference System:
User input: OSGB36 / British National Grid
wkt:
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]]st_crs(4326)$epsg[1] 4326st_crs(3857)Coordinate Reference System:
User input: EPSG:3857
wkt:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]st_crs(27700)$proj4string[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"cholera_4326 <- CholeraDeaths %>%
st_transform(4326)st_bbox(cholera_4326) xmin ymin xmax ymax
-0.1400738 51.5118557 -0.1329335 51.5158345 ggplot(cholera_4326) +
annotation_map_tile(type = "osm", zoomin = 0) +
geom_sf(aes(size = Count), alpha = 0.7)Zoom: 17
help("spTransform-methods", package = "sp")st_crs(CholeraDeaths)$proj4string[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"cholera_latlong <- CholeraDeaths %>%
st_set_crs(27700) %>%
st_transform(4326)
snow <- ggplot(cholera_latlong) +
annotation_map_tile(type = "osm", zoomin = 0) +
geom_sf(aes(size = Count))pumps <- st_read(dsn, layer = "Pumps")Reading layer `Pumps' from data source
`/home/esuess/classes/2023-2024/Fall 2023/Stat651/Presentations/08_spatial/spatial4/SnowGIS_SHP'
using driver `ESRI Shapefile'
Simple feature collection with 8 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 529183.7 ymin: 180660.5 xmax: 529748.9 ymax: 181193.7
Projected CRS: OSGB36 / British National Gridpumps_latlong <- pumps %>%
st_set_crs(27700) %>%
st_transform(4326)
snow +
geom_sf(data = pumps_latlong, size = 3, color = "red")Zoom: 17
Leaflet
library(tidygeocoder)
Attaching package: 'tidygeocoder'The following object is masked from 'package:prettymapr':
geocodewhite_house <- tibble(
address = "The White House, Washington, DC"
) %>%
tidygeocoder::geocode(address, method = "osm")Passing 1 address to the Nominatim single address geocoderQuery completed in: 1 secondslibrary(leaflet)
white_house_map <- leaflet() %>%
addTiles() %>%
addMarkers(data = white_house)Assuming "long" and "lat" are longitude and latitude, respectivelywhite_house_map white_house <- white_house %>%
mutate(
title = "The White House",
street_address = "1600 Pennsylvania Ave"
)
white_house_map %>%
addPopups(
data = white_house,
popup = ~paste0("<b>", title, "</b></br>", street_address)
)Assuming "long" and "lat" are longitude and latitude, respectively