Working with spatial data

Peter Ralph

4 March 2021 – Advanced Biological Statistics

An example

Biketown

git clone git@github.com:UO-Biostats/datasets.git

catmapper: Wikimedia

Steve Morgan: Wikimedia

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:

Starting and ending locations of trips.
Net flux of bikes by area (e.g., places where more trips start than end).
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.

library(sp)
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