Simple Features in R

Mark Cherrie
29/06/18

What are we doing today?

  • Introduction (Geospatial analysis in R, What are simple features?, Exposome and nudge theory)
  • Method (Setup, Google API's, Moves App, Secondary Natural Environment data)
  • Analysis (Activity Count, Hour Count, Minimal Convex Polygon, Bearing, Processing Google Streetview and Vision, Calculated Semantic Naturalness)
  • Conclusions (Share visualisations, Summary, next steps, references)
  • Have a go!

Geospatial Analysis in R

What are Simple Features?

  • Simple feature access is an ISO standard that is widely adopted. It is used in spatial databases (PostGIS), GIS (ArcGIS), open source libraries, GeoJSON, GeoSPARQL, etc.
    • Fast reading and writing of data
    • Enhanced plotting performance
    • sf objects can be treated as data frames in most operations
    • sf functions can be combined using %>% operator and works well with the tidyverse collection of R packages

Exposome and Nudge Theory

Exposome and Nudge Theory

Data Science Pipeline

Setup

  • Install external packages
library(devtools)
#install_github("r-spatial/sf")
#install.packages("RoogleVision", repos = c(getOption("repos"), "http://cloudyr.github.io/drat"))
#install_github("r-spatial/mapview@develop")
  • Install CRAN packages
#install.packages(c("tidyverse",
# "data.table", "pbapply","XML", 
# "dtplyr","sp", "googleway",
# "adehabitatHR", "zoo", "argosfilter"))

Setup

  • Load packages
invisible(lapply(c(
         "sf","tidyverse","data.table",
         "pbapply","XML", "dtplyr", 
         "mapview","sp", "googleway", 
         "RoogleVision", "adehabitatHR", 
         "zoo", "argosfilter"), 
         require, character.only = TRUE))
  • Load external functions
# Install User Written Functions
for (i in c("functions")){
  source(paste0(i, ".R"), echo=FALSE)
}

Google API's

  • There is lots of information available in Google application programming interfaces:

    • Places API
    • Directions API
    • Street View API
    • Vision API
  • Get an API key: https://developers.google.com/maps/documentation/javascript/get-api-key

  • Note that starting July 16, 2018, a new pay-as-you-go pricing plan will go into effect for Maps, Routes, and Places ($200 free a month; e.g. 40,000 directions calls)

Google Maps Data

  • Get route polyline
Usertrack<-GoogleRoute(mapkey= "AIzaSyA8d6ToR2dww8BcfrYhiefkhB0DoZZU8Ck", latO=55.934544, lonO=-3.228447, latD=55.947779, lonD=-3.184145, mode="walking")
plot(Usertrack)

plot of chunk unnamed-chunk-2

Moves App Data

tracks<-get_tracks("data/moves_export/")

activi<-get_activities("data/moves_export/")

places<-get_places("data/moves_export/")
  • Note we'll concentrate on data from tracks today

Natural Environment Data

Natural Environment Data

  • sf::st_read() is substantially faster than rgdal::readOGR()
    • Approx. 14-18 times faster for Shapefile and KML files
library(sf)
path="data/opgrsp_essh_nt/OS Open Greenspace (ESRI Shape File) NT/data/"

greenspaces <- read_sf(paste0(path,"NT_GreenspaceSite.shp"))

greenaccesspoints <- 
read_sf(paste0(path, "NT_AccessPoint.shp"))

quietwalks <- read_sf("data/Quiet_routes/Quiet_routes.kml")
  • Note read_sf instead of st_read

Plotting Simple Features

  • Simple plot
tracks %>%
  st_as_sf(coords = c("longitude","latitude")) %>%
  st_set_crs(4326) %>%
  st_cast("MULTIPOINT") %>%
  plot(.)
  • Note sf function names are relatively consistent and intuitive (all begin with st_)
  • Note that you can set the coordinate reference system using EPSG codes (http://epsg.io)

Plotting Simple Features

plot of chunk unnamed-chunk-6

Plotting Simple Features - Mapview

# Options
mapviewOptions(basemaps = c("CartoDB.Positron", "OpenStreetMap","Esri.WorldImagery"), layers.control.pos = "topright")

# Plot
m1  <-  mapview(quietwalks,  zcol = "Name")  +
        mapview(greenspaces, zcol = "function.", alpha = 0) + 
        mapview(greenaccesspoints, zcol="accessType", alpha = 0)

# Static Output
# mapshot(m1, file = "quietwalks.png")
  • Note that alpha=0 takes out the polygon boundary
  • Note that the mapshot function is useful if you need a static image

Plotting Simple Features - Mapview

Pipes - Activity Count

# Atribute Data- Activity count
tracks %>%
  st_as_sf(coords = c("longitude","latitude")) %>%
  st_set_crs(4326) %>%
  st_cast("MULTIPOINT") %>%
  mutate(count=1) %>%
  group_by(activity) %>%
  summarise(activitycount=sum(count)) %>%
  ggplot(., aes(x=activity, y=activitycount)) + geom_bar(stat = "identity")
  • Note that %>% will pipe a value forward into an expression or function call
  • Use . if you want to take the manipulated dataframe

Pipes - Activity Count

plot of chunk unnamed-chunk-9

Pipes - Hour Count

# Time data - Hour count
tracks %>%
  sf::st_as_sf(coords = c("longitude","latitude")) %>%
  st_set_crs(4326) %>%
  st_cast("MULTIPOINT") %>%
  mutate(count=1) %>%
  mutate(hour=format(as.POSIXct(strptime(time,"%Y-%m-%dT%H:   %M",tz="GMT")),format = "%H")) %>%
  group_by(hour) %>%
  summarise(hourcount=sum(count)) %>%
  ggplot(., aes(x=hour, y=hourcount)) + geom_bar(stat = "identity")
  • If you think you will have loaded packaged using the same functions, then use 'packagename::' to make it explicit that the function is from a specific package

Pipes - Hour Count

plot of chunk unnamed-chunk-11

Pipes - Minimal Convex Polygon

# Geographic data - calculate Minimal convex polygon for each 
tracks_sf_pt<-tracks %>%
  sf::st_as_sf(coords = c("longitude","latitude")) %>%
  sf::st_set_crs(4326) %>%
  sf::st_transform(., 27700)
tracksSPDF<-as(tracks_sf_pt, "Spatial")
track_poly <- mcp(tracksSPDF[,4], percent=95)
mcp_sf_pt<-st_as_sf(track_poly, coords = c("longitude", "latitude"),crs = 27700)
m2<-mapview(mcp_sf_pt, zcol = "area", alpha.regions = 0.3)
  • st_distance(), st_length() and st_area() will report results with values in the current CRS, so we convert to British National Grid projection for output in metres.
  • Conversion of simple feature objects of class sf or sfc into corresponding Spatial* objects is done using the as method, coercing to Spatial

Pipes - Minimal Convex Polygon

Pipes - Bearing

# Calculate bearing
# convert to spatialMultipointsdataframe
tracks_sf_pt_bearing<-tracks_sf_pt %>%
  st_cast("MULTIPOINT") %>%
  filter(date=="2016-09-01") %>%
  as(., "Spatial")
tracks_sf_pt_bearing<-as.data.frame(tracks_sf_pt_bearing)
bearing<-bearingTrack(tracks_sf_pt_bearing$X1, tracks_sf_pt_bearing$X2)
bearing<-append(bearing, 90) 
tracks_sf_pt_bearing<-cbind(tracks_sf_pt_bearing, bearing)
  • Note that is also possible to create data.frame objects with geometry list-columns that are not of class sf

Fill Missing and Format

# fill in missing from last unmissing
tracks_sf_pt_bearing$bearing<-zoo::na.locf(tracks_sf_pt_bearing$bearing)
# Specify degrees clockwise from true north, around the camera locus
tracks_sf_pt_bearing$bearing<-ifelse(tracks_sf_pt_bearing$bearing<0, tracks_sf_pt_bearing$bearing+360, tracks_sf_pt_bearing$bearing+0)
# Format
colnames(tracks_sf_pt_bearing)[1:2]<-c("longitude", "latitude")

# Add bearings to 45 degrees to the left and right of the image
tracks_sf_pt_bearingP45<-tracks_sf_pt_bearing
tracks_sf_pt_bearingM45<-tracks_sf_pt_bearing
tracks_sf_pt_bearingP45$bearing<-tracks_sf_pt_bearing$bearing+45
tracks_sf_pt_bearingM45$bearing<-tracks_sf_pt_bearing$bearing-45
tracks_sf_pt_bearing<-rbind(tracks_sf_pt_bearing, tracks_sf_pt_bearingP45,tracks_sf_pt_bearingM45)

Streetview Batch Processing

# Create a new directory for the streetview images if one doesn't already exist
dir.create(file.path(getwd(), "streetview_images"), showWarnings = FALSE)

# Streetview function
imagedownloader<-function(latitude,longitude, bearing){
  png(paste0("streetview_images/",latitude,"_", longitude,"_", bearing, ".png"), width=640, height=480)
  google_streetview(location = c(latitude,longitude), size = c(640,480), panorama_id = NULL, output = "plot", heading = bearing, fov = 90, pitch = 0, response_check = FALSE, key = "AIzaSyA8d6ToR2dww8BcfrYhiefkhB0DoZZU8Ck")
  dev.off()
}
  • R can perform system commands (e.g. create directory)

Streetview Batch Processing

# Batch download images
tracks_sf_pt_bearing<-subset(tracks_sf_pt_bearing, select=c("latitude", "longitude", "bearing"))
library(plyr)
mdply(tracks_sf_pt_bearing, imagedownloader)
  • One of the most useful functions is mdply - a multi-argument function with values taken from columns of an data frame or array, which combines results into a data frame

Streetview Batch Processing

  • Note that we need to manually check the pictures, some of them are inside, some of them have things in the way, for now we'll just need to manaully delete.