Spatial Data

Spatial data

library(tidyverse)
library(mdsr)
library(sp)
plot(CholeraDeaths)

It is probably best to do this in an R Project.

library(rgdal)

rgdal: version: 1.3-6, (SVN revision 773)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
 Path to GDAL shared files: /usr/share/gdal/2.1
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.3-1 

Attaching package: ‘rgdal’

The following object is masked from ‘package:mosaic’:

    project

download.file("http://rtwilson.com/downloads/SnowGIS_SHP.zip",
              dest="SnowGIS.zip", mode="wb")

trying URL 'http://rtwilson.com/downloads/SnowGIS_SHP.zip'
Content type 'application/zip' length 6875824 bytes (6.6 MB)
==================================================
downloaded 6.6 MB

unzip("SnowGIS.zip")
getwd()

[1] "/home/esuess/classes/2018-2019/01 - Fall 2018/Stat651/Presentations/08_spatial"

dsn <- paste0("./SnowGIS_SHP/")
list.files(dsn)

 [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"            

ogrListLayers(dsn)

[1] "Pumps"          "Cholera_Deaths"
attr(,"driver")
[1] "ESRI Shapefile"
attr(,"nlayers")
[1] 2

ogrInfo(dsn, layer = "Cholera_Deaths")

Source: "/home/esuess/classes/2018-2019/01 - Fall 2018/Stat651/Presentations/08_spatial/SnowGIS_SHP", layer: "Cholera_Deaths"
Driver: ESRI Shapefile; number of rows: 250 
Feature type: wkbPoint with 2 dimensions
Extent: (529160.3 180857.9) - (529655.9 181306.2)
CRS: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs  
LDID: 87 
Number of fields: 2 
   name type length typeName
1    Id    0      6  Integer
2 Count    0      4  Integer

CholeraDeaths <- readOGR(dsn, layer = "Cholera_Deaths")

OGR data source with driver: ESRI Shapefile 
Source: "/home/esuess/classes/2018-2019/01 - Fall 2018/Stat651/Presentations/08_spatial/SnowGIS_SHP", layer: "Cholera_Deaths"
with 250 features
It has 2 fields

summary(CholeraDeaths)

Object of class SpatialPointsDataFrame
Coordinates:
               min      max
coords.x1 529160.3 529655.9
coords.x2 180857.9 181306.2
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs]
Number of points: 250
Data attributes:
       Id        Count       
 Min.   :0   Min.   : 1.000  
 1st Qu.:0   1st Qu.: 1.000  
 Median :0   Median : 1.000  
 Mean   :0   Mean   : 1.956  
 3rd Qu.:0   3rd Qu.: 2.000  
 Max.   :0   Max.   :15.000  

str(CholeraDeaths@data)

'data.frame':   250 obs. of  2 variables:
 $ Id   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Count: int  3 2 1 1 4 2 2 2 3 2 ...

Making maps

cholera_coords <- as.data.frame(coordinates(CholeraDeaths))
cholera_coords

cholera_coords %>% ggplot(aes(x = coords.x1, y = coords.x2)) +
  geom_point() +
  coord_quickmap()

library(ggmap)
# I have saved my GOOGLEKEY as a system variable.  
# If you are interested in this you can ask about it in office hours.
register_google(key = Sys.getenv("GOOGLEKEY"))
# register_google(key = "AIzaSyCtAkBsI2ZlnrxVyTiiOrp")  # This is the line to add your 
                                                        # Google API key and uncomment to run.
# google_key()  # The uncomment this line to see if your key has been read into R ok.
m <- get_map("John Snow, London, England", zoom = 17, maptype = "roadmap") 

Source : https://maps.googleapis.com/maps/api/staticmap?center=John%20Snow,%20London,%20England&zoom=17&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-C8Qslw
Source : https://maps.googleapis.com/maps/api/geocode/json?address=John%20Snow%2C%20London%2C%20England&key=xxx-C8Qslw

ggmap(m)

# Alternatively use latlong.net to get the coordinates.
m <- get_map(location = c(lon = -0.136530, lat = 51.513378), zoom = 17, maptype = "roadmap")

Source : https://maps.googleapis.com/maps/api/staticmap?center=51.513378,-0.13653&zoom=17&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-C8Qslw

ggmap(m)

As the book says, you will see nothing on the map.

ggmap(m) + geom_point(data = as.data.frame(CholeraDeaths),
                      aes(x = coords.x1, y = coords.x2, size= Count))

head(as.data.frame(CholeraDeaths))

Using different units.

attr(m, "bb")

Projections

library(maps)

Attaching package: ‘maps’

The following object is masked from ‘package:purrr’:

    map

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)

proj4string(CholeraDeaths) %>% strwrap()

[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"

European Petroleum Survery Group (EPSG) code provide a shorthand for the full PROJ.4 strings.

Coordinate Reference System (CRS)

CRS("+init=epsg:4326")

CRS arguments:
 +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

CRS("+init=epsg:3857")

CRS arguments:
 +init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m
+nadgrids=@null +no_defs 

CRS("+init=epsg:27700")

CRS arguments:
 +init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36
+units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 

cholera_latlong <- CholeraDeaths %>% spTransform(CRS("+init=epsg:4326"))

Now looks to be longitude and latatude.

head(cholera_latlong)

             coordinates Id Count
1 (-0.1363245, 51.51291)  0     3
2 (-0.1362775, 51.51285)  0     2
3 (-0.1362473, 51.51281)  0     1
4 (-0.1362063, 51.51275)  0     1
5 (-0.1361612, 51.51269)  0     4
6 (-0.1359313, 51.51267)  0     2

Plot the map, but note the points do not look right. Some point are in the street.

ggmap(m) + geom_point(aes(x = coords.x1, y = coords.x2, size = Count), 
                      data = as.data.frame(cholera_latlong))

The authors examine the meta data and fix the problem

help("spTransform-methods", package = "rgdal")

CholeraDeaths %>% proj4string() %>% showEPSG()

[1] "OGRERR_UNSUPPORTED_SRS"

proj4string(CholeraDeaths) <- CRS("+init=epsg:27700")

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=m +no_defs
without reprojecting.
For reprojection, use function spTransform

cholera_latlong <- CholeraDeaths %>% 
  spTransform(CRS("+init=epsg:4326")) 
snow <- ggmap(m) + geom_point(aes(x = coords.x1, y = coords.x2,
                                  size = Count), data = as.data.frame(cholera_latlong))
pumps <- readOGR(dsn, layer = "Pumps")

OGR data source with driver: ESRI Shapefile 
Source: "/home/esuess/classes/2018-2019/01 - Fall 2018/Stat651/Presentations/08_spatial/SnowGIS_SHP", layer: "Pumps"
with 8 features
It has 1 fields

proj4string(pumps) <- CRS("+init=epsg:27700") 

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=m +no_defs
without reprojecting.
For reprojection, use function spTransform

pumps_latlong <- pumps %>% 
  spTransform(CRS("+init=epsg:4326")) 
snow + geom_point(data = as.data.frame(pumps_latlong),
                  aes(x = coords.x1, y = coords.x2), size = 3, color = "red")

Geocoding, routs, and distances

smith <- "Smith College, Northampton, MA 01063" 
geocode(smith)

Source : https://maps.googleapis.com/maps/api/geocode/json?address=Smith%20College%2C%20Northampton%2C%20MA%2001063&key=xxx-C8Qslw

amherst <- "Amherst College, Amherst, MA"
geocode(amherst)

Source : https://maps.googleapis.com/maps/api/geocode/json?address=Amherst%20College%2C%20Amherst%2C%20MA&key=xxx-C8Qslw

mapdist(from = smith, to = amherst, mode = "driving")

Source : https://maps.googleapis.com/maps/api/distancematrix/json?origins=Smith%20College%2C%20Northampton%2C%20MA%2001063&destinations=Amherst%20College%2C%20Amherst%2C%20MA&mode=driving&language=en-EN&key=xxx-C8Qslw
REQUEST_DENIED : This API project is not authorized to use this API.matching was not perfect, returning what was found.
Error in `*tmp*`[[c(1, 1)]] : no such index at level 1

legs_df <- route(smith, amherst, alternatives = TRUE) 

Source : https://maps.googleapis.com/maps/api/directions/json?origin=Smith%20College%2C%20Northampton%2C%20MA%2001063&destination=Amherst%20College%2C%20Amherst%2C%20MA&mode=driving&units=metric&alternatives=true&key=xxx-C8Qslw
Error in rep(LETTERS[k], stepsPerRoute[k]) : invalid 'times' argument

qmap("The Quarters, Hadley, MA", zoom = 12, maptype = 'roadmap') + 
  geom_leg(aes(x = startLon, y = startLat, xend = endLon, yend = endLat), 
           alpha = 3/4, size = 2, color = "blue", data = legs_df)

Source : https://maps.googleapis.com/maps/api/staticmap?center=The%20Quarters,%20Hadley,%20MA&zoom=12&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-C8Qslw
Source : https://maps.googleapis.com/maps/api/geocode/json?address=The%20Quarters%2C%20Hadley%2C%20MA&key=xxx-C8Qslw
Error in fortify(data) : object 'legs_df' not found

Leaflet

white_house <- geocode("The White House, Washington, DC") 

Source : https://maps.googleapis.com/maps/api/geocode/json?address=The%20White%20House%2C%20Washington%2C%20DC&key=xxx-C8Qslw

white_house

library(leaflet) 
map <- leaflet() %>%
  addTiles() %>% 
  addMarkers(lng = ~lon, lat = ~lat, data = white_house)
white_house <- white_house %>% mutate(title = "The White House", 
                  address = "2600 Pennsylvania Ave")
map %>% addPopups(lng = ~lon, lat = ~lat, data = white_house, 
                  popup = ~paste0("<b>", title, "</b></br>", address))

RDSTK

An alternative package for geocoding.

library(RDSTK)

Loading required package: plyr
---------------------------------------------------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
---------------------------------------------------------------------------------------------------------------------------

Attaching package: ‘plyr’

The following object is masked from ‘package:maps’:

    ozone

The following object is masked from ‘package:mosaic’:

    count

The following objects are masked from ‘package:plotly’:

    arrange, mutate, rename, summarise

The following objects are masked from ‘package:dplyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

The following object is masked from ‘package:purrr’:

    compact

Loading required package: rjson
Loading required package: RCurl
Loading required package: bitops

Attaching package: ‘RCurl’

The following object is masked from ‘package:tidyr’:

    complete

# alternative to geocode from ggmap package for geocoding using OpenStreeMap
white_house <- street2coordinates("1600 Pennsylvania Avenue, Washington, DC")
white_house

map <- leaflet() %>%
  addTiles() %>% 
  addMarkers(lng = ~longitude, lat = ~latitude, data = white_house)
white_house <- white_house %>% mutate(title = "The White House", 
                  address = "2600 Pennsylvania Ave")
map %>% addPopups(lng = ~longitude, lat = ~latitude, data = white_house, 
                  popup = ~paste0("<b>", title, "</b></br>", address))

Extended example: Historical airline route maps

library(nycflights13)
my_carrier <- "DL"
my_year <- 2013
airlines

airports

flights

destinations <- flights %>% filter(year == my_year, carrier == my_carrier) %>% 
  left_join(airports, by = c("dest" = "faa")) %>% 
  group_by(dest) %>%
  summarise( N=n(), lon = max(lon), lat = mean(lat) ) %>%
  collect()
glimpse(destinations)

Observations: 40
Variables: 4
$ dest <chr> "ATL", "AUS", "BNA", "BOS", "BUF", "CVG", "DCA", "DEN", "DTW", "EYW", "FLL", "IND", "JAC", "JAX", "LAS", "L...
$ N    <int> 10571, 357, 1, 972, 3, 4, 2, 1043, 3875, 17, 2903, 2, 2, 1, 1673, 2501, 82, 3663, 432, 2929, 2864, 1129, 1,...
$ lon  <dbl> -84.42807, -97.66989, -86.67819, -71.00518, -78.73217, -84.66782, -77.03772, -104.67318, -83.35339, -81.759...
$ lat  <dbl> 33.63672, 30.19453, 36.12447, 42.36435, 42.94053, 39.04884, 38.85208, 39.86166, 42.21244, 24.55611, 26.0725...

segments <- flights %>% filter(year == my_year, carrier == my_carrier) %>% 
  group_by(origin, dest) %>% summarize(N = n()) %>% 
  left_join(airports, by = c("origin" = "faa")) %>% 
  left_join(airports, by = c("dest" = "faa")) %>% collect() %>% 
  na.omit()
dim(segments)

[1] 53 17

route_map <- qmap("junction city, kansas", zoom = 4, maptype = "roadmap") +
  geom_point(data = destinations, alpha = 0.5, aes(x = lon, y = lat, size = N)) +
  scale_size() +
  theme_map() 

Source : https://maps.googleapis.com/maps/api/staticmap?center=junction%20city,%20kansas&zoom=4&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-C8Qslw
Source : https://maps.googleapis.com/maps/api/geocode/json?address=junction%20city%2C%20kansas&key=xxx-C8Qslw

route_map

destinations %>% arrange(desc(N))

route_map + geom_segment( 
  aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, color = N), 
  size = 0.05, arrow = arrow(length = unit(0.3, "cm")), data = segments)

lines <- bind_rows( segments %>%
  select(origin, dest, lat.x, lon.x) %>%
  rename(lat = lat.x, lon = lon.x), segments %>%
  select(origin, dest, lat.y, lon.y) %>%
  rename(lat = lat.y, lon = lon.y)) %>% 
  arrange(origin, dest) %>% na.omit() %>% 
  group_by(origin, dest) %>%
  do(line = Line(as.data.frame(select(., lon, lat))) 
)    
head(lines, 3)

make_line <- function(x) { 
  Lines(list(x[["line"]]), ID = paste0(x$origin, "-", x$dest))
}
lines_list <- apply(lines, MARGIN = 1, make_line)

segments_sp <- SpatialLines(lines_list, CRS("+proj=longlat")) 
summary(segments_sp)

Object of class SpatialLines
Coordinates:
         min       max
x -122.59750 -70.30928
y   24.55611  47.44900
Is projected: FALSE 
proj4string : [+proj=longlat +ellps=WGS84]

segments_sp <- segments_sp %>% spTransform(CRS("+init=epsg:4326"))

library(leaflet) 
l_map <- leaflet() %>%
  addTiles() %>% 
  addCircles(lng = ~lon, lat = ~lat, weight = 1, radius = ~sqrt(N) * 500, data = destinations) %>% 
  addPolylines(weight = 0.4, data = segments_sp) %>% setView(lng = -80, lat = 38, zoom = 6)

Data contains 2 rows with either missing or invalid lat/lon values and will be ignored

l_map