\[ %% % Add your macros here; they'll be included in pdf and html output. %% \newcommand{\R}{\mathbb{R}} % reals \newcommand{\E}{\mathbb{E}} % expectation \renewcommand{\P}{\mathbb{P}} % probability \DeclareMathOperator{\logit}{logit} \DeclareMathOperator{\logistic}{logistic} \DeclareMathOperator{\SE}{SE} \DeclareMathOperator{\sd}{sd} \DeclareMathOperator{\var}{var} \DeclareMathOperator{\cov}{cov} \DeclareMathOperator{\cor}{cor} \DeclareMathOperator{\Normal}{Normal} \DeclareMathOperator{\MVN}{MVN} \DeclareMathOperator{\LogNormal}{logNormal} \DeclareMathOperator{\Poisson}{Poisson} \DeclareMathOperator{\Beta}{Beta} \DeclareMathOperator{\Binom}{Binomial} \DeclareMathOperator{\Gam}{Gamma} \DeclareMathOperator{\Exp}{Exponential} \DeclareMathOperator{\Cauchy}{Cauchy} \DeclareMathOperator{\Unif}{Unif} \DeclareMathOperator{\Dirichlet}{Dirichlet} \DeclareMathOperator{\Wishart}{Wishart} \DeclareMathOperator{\StudentsT}{StudentsT} \DeclareMathOperator{\Weibull}{Weibull} \newcommand{\given}{\;\vert\;} \]

Working with spatial data

Peter Ralph

Advanced Biological Statistics

An example

Biketown

git clone git@github.com:UO-Biostats/datasets.git
Albert in a bike (catmapper:Wikimedia)

catmapper: Wikimedia

biketown bikes (Steve Morgan: Wikimedia)

Steve Morgan: Wikimedia

how biketown works: screenshot

August, 2018

data_dir <- "../../../datasets" 
bt <- read.csv(file.path(data_dir, "biketown/2019_08.csv.gz"), header=TRUE, stringsAsFactors=FALSE)
xylims <- c(-122.75, -122.6, 45.46, 45.56)
bt <- subset(bt, !is.na(Start_Latitude) & Start_Latitude > xylims[3] & Start_Latitude < xylims[4] 
                 & !is.na(End_Latitude) & End_Latitude > xylims[3] & End_Latitude < xylims[4] 
                 & !is.na(Start_Longitude) & Start_Longitude > xylims[1] & Start_Longitude < xylims[2] 
                 & !is.na(End_Longitude) & End_Longitude > xylims[1] & End_Longitude < xylims[2])
nrow(bt)
## [1] 38587
layout(t(1:2))
plot(Start_Latitude ~ Start_Longitude, data=bt, pch=20, cex=0.5, col=adjustcolor('black', 0.2), xlab='eastings', ylab='northings')
plot(Start_Latitude ~ Start_Longitude, data=bt, type='n', xlab='eastings', ylab='northings')
with(bt, segments(x0=Start_Longitude, y0=Start_Latitude,
                  x1=End_Longitude, y1=End_Latitude, col=adjustcolor("black", 0.05), lwd=0.5))

plot of chunk r biketown_it

Goals:

Visualize:

  1. Starting and ending locations of trips.
  2. Net flux of bikes by area (e.g., places where more trips start than end).
  3. Distance of trips depending on start and end location.

Making nice maps

See this tutorial by Manual Gimond.

library(rgdal)
library(sp)
library(raster)

How I got the mapping data

crop_extent <- extent(xylims)
river <- readOGR(dsn=file.path(data_dir, "biketown/River_Overlay/"))
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/peter/teaching/uo-biostats/datasets/biketown/River_Overlay", layer: "River_Overlay_Boundary"
## with 7 features
## It has 6 fields
bikepaths <- readOGR(dsn=file.path(data_dir, "biketown/bicycle_network/"))
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/peter/teaching/uo-biostats/datasets/biketown/bicycle_network", layer: "Bicycle_Network"
## with 39798 features
## It has 10 fields
streets <- readOGR(dsn=file.path(data_dir, "biketown/pdx_streets/"))
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/peter/teaching/uo-biostats/datasets/biketown/pdx_streets", layer: "Streets"
## with 108454 features
## It has 39 fields

How I got the mapping data, continued

bikepaths <- subset(bikepaths, Status == "Active")
streets <- subset(streets, LCITY == "Portland" | RCITY == "Portland")
river <- crop(river, crop_extent)
freeways <- crop(subset(streets, TYPE == 1110), crop_extent)
bigstreets <- crop(subset(streets, TYPE %in% c(1300, 1400)), crop_extent)
rm(streets)

How you can get the mapping data

show(load(file.path(data_dir, "biketown/pdx_features.RData"))) 
## [1] "crop_extent" "river"       "bikepaths"   "freeways"    "bigstreets" 

layout(t(1:2))
sp::plot(river, col=adjustcolor('blue', 0.25), main='river')
sp::plot(bikepaths, col=adjustcolor('grey', 0.75), main='bikepaths')

plot of chunk r plot_river

bt_start <- SpatialPointsDataFrame(coords=cbind(bt$Start_Longitude, bt$Start_Latitude),
                                   data=bt, proj4string=CRS("+proj=longlat"))
bt_end <- SpatialPointsDataFrame(coords=cbind(bt$End_Longitude, bt$End_Latitude),
                                   data=bt, proj4string=CRS("+proj=longlat"))
bt_start <- crop(
                 spTransform(bt_start,
                             proj4string(river)),
                 crop_extent)
bt_end <- crop(
               spTransform(bt_end,
                           proj4string(river)),
               crop_extent)

Important: the “projection” tells R what the coordinates mean:

bt_start <- SpatialPointsDataFrame( # make a legit 'spatial' object
    coords=cbind(bt$Start_Longitude, bt$Start_Latitude), # at these coordinates
    data=bt,
    proj4string=CRS("+proj=longlat") # <---- the projection!
)

… and the proj4 string is the magical incantation that says what the projection is.

setup_map <- function (..., add=FALSE) {
    plot(bigstreets, col=grey(0.75), add=add, ...)
    lines(freeways, col=adjustcolor('brown', 0.25), lwd=2)
    plot(river, col=adjustcolor('blue', 0.25), add=TRUE)
    lines(bikepaths, col=adjustcolor(grey(0.25), 0.5))
}

layout(t(1:2))
setup_map(main='starts'); points(bt_start, pch=20, cex=0.5, col='red')
setup_map(main='ends'); points(bt_end, pch=20, cex=0.5, col='blue')

plot of chunk r plot_pdx

More resources:

// reveal.js plugins