We can import the envirocar directly into spatial objects by exploiting the GeoJSON format which is supported by OGR, hence R function readOGR in package rgdal. To read from a https connection, one needs RCurl when on windows.
url = "https://giv-car.uni-muenster.de/stable/rest/tracks/5207d871e4b058cd3d669afe"
require(rgdal) # readOGR
## Loading required package: rgdal Loading required package: sp rgdal:
## version: 0.8-11, (SVN revision 479M) Geospatial Data Abstraction Library
## extensions to R successfully loaded Loaded GDAL runtime: GDAL 1.10.0,
## released 2013/04/24 Path to GDAL shared files: /usr/share/gdal/1.10 Loaded
## PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] Path to PROJ.4
## shared files: (autodetected)
layer <- readOGR(url, layer = "OGRGeoJSON")
## OGR data source with driver: GeoJSON
## Source: "https://giv-car.uni-muenster.de/stable/rest/tracks/5207d871e4b058cd3d669afe", layer: "OGRGeoJSON"
## with 173 features and 3 fields
## Feature type: wkbPoint with 2 dimensions
class(layer)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
layer[1, ]
## coordinates id time
## 1 (7.59584, 51.9691) 0 2013-08-11T16:56:06Z
## phenomenons
## 1 { "Intake Temperature": { "value": 38.0, "unit": "c" }, "Intake Pressure": { "value": 126.0, "unit": "kPa" }, "CO2": { "value": 0.00587494280504, "unit": "kg\\/h" }, "Calculated MAF": { "value": 26.7454, "unit": "g\\/s" }, "Speed": { "value": 4.0, "unit": "km\\/h" }, "Rpm": { "value": 1339.0, "unit": "u\\/min" }, "Consumption": { "value": 0.002216959549072, "unit": "l\\/h" } }
As we can see, the attribute contains id (point nr) and time (still as factor, not as POSIXt) and the measured values are still in a JSON string that needs to be parsed. For this, we wrote an import function that converts time, parses the JSON into measured values, and adds a units attribute to the sp object:
importEnviroCar = function(file) {
require(rjson) # fromJSON
require(maptools) # spCbind
# read data as spatial object:
layer = readOGR(file, layer = "OGRGeoJSON")
# convert time from text to POSIXct:
layer$time = as.POSIXct(layer$time, format = "%Y-%m-%dT%H:%M:%SZ")
# the third column is JSON, we want it in a table (data.frame) form: 1. form
# a list of lists
l1 = lapply(as.character(layer[[3]]), fromJSON)
# 2. parse the $value elements in the sublist:
l2 = lapply(l1, function(x) as.data.frame(lapply(x, function(X) X$value)))
# bind these value elements, row by row, together
ret = do.call(rbind, l2)
# read the units:
units = lapply(l1[1], function(x) as.data.frame(lapply(x, function(X) X$unit)))
# add a units attribute to layer
layer[[3]] = NULL
# add the table as attributes to the spatial object
if (length(layer) == nrow(ret)) {
layer = spCbind(layer, ret)
attr(layer, "units") = units[[1]]
layer
} else NULL
}
url = "https://giv-car.uni-muenster.de/stable/rest/tracks/5207d871e4b058cd3d669afe"
sp = importEnviroCar(url)
## Loading required package: rjson Loading required package: maptools
## Checking rgeos availability: TRUE
## OGR data source with driver: GeoJSON
## Source: "https://giv-car.uni-muenster.de/stable/rest/tracks/5207d871e4b058cd3d669afe", layer: "OGRGeoJSON"
## with 173 features and 3 fields
## Feature type: wkbPoint with 2 dimensions
sp[1, ]
## coordinates id time Intake.Temperature
## 1 (7.59584, 51.9691) 0 2013-08-11 16:56:06 38
## Intake.Pressure CO2 Calculated.MAF Speed Rpm Consumption
## 1 126 0.005875 26.75 4 1339 0.002217
attr(sp, "units")
## Intake.Temperature Intake.Pressure CO2 Calculated.MAF Speed Rpm
## 1 c kPa kg/h g/s km/h u/min
## Consumption
## 1 l/h
library(lattice)
trellis.par.set(sp.theme())
spplot(sp, "Consumption", colorkey = TRUE)
Here, we interpolate the point measurements to a new set of points, using inverse distance weighting
n = length(sp) * 5 # create 5 points between each observation...
cc = coordinates(sp)
newx = approx(cc[, 1], n = n)$y # along x,
newy = approx(cc[, 2], n = n)$y # along y,
crs = CRS(proj4string(sp))
newpts = SpatialPoints(cbind(newx, newy), crs)
# Alternatively: convert to SpatialLines, then use spsample
library(gstat)
# interpolate Consumption values:
idw = idw(Consumption ~ 1, sp, newpts)
## [inverse distance weighted interpolation]
We can plot this sequence by
spplot(idw[1], colorkey = TRUE, scales = list(draw = TRUE), main = "inverse distance interpolated Fuel Consumption values")
Here, the new points are evenly distributed between all existing points, meaning that with 5 times as much points, between every two points 4 points are added equidistant. This gives a new point density that resembles the original point density.
Alternatively, we could generate interpolation points
We can aggregate to a 5 x 5 grid, using an arbitray aggregation function, here max
bb = bbox(sp)
grd = GridTopology(bb[, 1], apply(bb, 1, diff)/5, c(5, 5))
sp.agg = aggregate(sp[-1], SpatialGrid(grd, crs), max)
and show the results, for
spplot(sp.agg[-c(1, 2, 8)], colorkey = TRUE, sp.layout = list("sp.points", sp,
col = 3), main = "maximum measured value per grid cell")
When looking at temporal variability, we see that for this trajectory the values are very unequally distrubed, indicating a break in the trajectory
plot(y = sp$Consumption, x = sp$time)
title("Fuel consumption over time")
We can also aggregate to time intervals, e.g. 10 min values. For this, we will convert the points to a zoo object (time series)
library(zoo)
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
xx = zoo(sp$Consumption, sp$time)
plot(aggregate(xx, cut(sp$time, "5 min"), mean))
points(xx)
title("10-minute mean Consumption (l/h)")
After having played with spatial objects (package sp) and temporal objects (package zoo), we can convert these data into spatio-temporal objects (package spacetime):
library(spacetime)
stidf = STIDF(geometry(sp), sp$time, sp@data)
stplot(geometry(stidf), main = "trajectory points, cut in 6 time slices with equal amount of points")
plots the geometry of the spatial points, over time, cutting time in six parts with an equal number of observations.
Trajectory data can be represented as Track object,
track = Track(stidf)
tracks = Tracks(list(tr1 = track))
tracksCollection = TracksCollection(list(tr = tracks))
stplot(tracksCollection)
The next plot adds attribute values as colour:
stplot(tracksCollection, attr = "speed", lwd = 3, scales = list(draw = T))
Using package ggmap, we can add a google map background (or other background):
bb = matrix(NA, 2, 2)
bb[2, ] = c(51.94, 51.985)
bb[1, ] = c(7.58, 7.67)
library(ggmap)
## Loading required package: ggplot2
map4ao <- get_map(location = as.vector(bb))
## Warning: bounding box given to google - spatial extent only approximate.
## converting bounding box to center/zoom specification. (experimental) Map
## from URL :
## http://maps.googleapis.com/maps/api/staticmap?center=51.9625,7.625&zoom=13&size=%20640x640&scale=%202&maptype=terrain&sensor=false
## Google Maps API Terms of Service : http://developers.google.com/maps/terms
# maptype = 'satellite', scale=2, zoom=4) Read the attributes to pass them
# to grid.raster
bbMap <- attr(map4ao, "bb")
latCenter <- with(bbMap, ll.lat + ur.lat)/2
lonCenter <- with(bbMap, ll.lon + ur.lon)/2
height <- with(bbMap, ur.lat - ll.lat)
width <- with(bbMap, ur.lon - ll.lon)
## Another component of the sp.layout in spplot
sp.raster <- list("grid.raster", map4ao, x = lonCenter, y = latCenter, width = width,
height = height, default.units = "native")
stplot(tracksCollection, scales = list(draw = TRUE), sp.layout = sp.raster,
col = "red", lwd = 2, main = "google map background")
The following code allows conversion from SpatialPoints into SpatialLines, and has been added to sp:
setAs("SpatialPoints", "Line", function(from) Line(coordinates(from)))
setAs("SpatialPoints", "Lines", function(from) Lines(as(from, "Line"), "ID"))
setAs("SpatialPoints", "SpatialLines", function(from) SpatialLines(list(as(from,
"Lines")), from@proj4string))
We can generalize a trajectory using Douglas-Peuker by using gSimplify in rgeos:
sl = as(sp, "SpatialLines")
library(rgeos)
## rgeos version: 0.2-20, (SVN revision (unknown)) GEOS runtime version:
## 3.3.8-CAPI-1.7.8 Polygon checking: TRUE
plot(gSimplify(sl, 5e-04), axes = TRUE) # WRONG: rgeos takes lat/long as Euclidian coords
plot(sp, add = T, col = "red")
title("generalization in Long/Lat (wrong!)")
however, without warning this falsely assumes that coordinates are metric, i.e. in a Euclidian system. They are not:
proj4string(sp)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
How can we resolve this?
Package rgdal contains code to reproject data:
library(rgdal)
utm = CRS("+proj=utm +zone=32 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
slT = spTransform(sl, utm)
5e-04 * 110000 # approximate meters per degree latitude
## [1] 55
plot(gSimplify(slT, 55), axes = TRUE) # RIGHT: rgeos takes coords Euclidian
plot(spTransform(sp, utm), add = T, col = "red")
title("generalization in UTM zone 32N")
A plot of the (projected) simplified trajectory, here represented by a set of selected points, is created by 1. conversion to lines, 2. simplifying the line, 3. conversion to points, 4. matching to the full set of points, and 5. selection:
sl.simplified = gSimplify(slT, 55) # 2
sp.simplified = as(sl.simplified, "SpatialPoints") # 3
sp = spTransform(sp, utm)
sel = zerodist2(sp.simplified, sp, 0.001)[, 2] # 4
sp.selected = sp[sel, ] # 5
length(sp.selected)
## [1] 16
sp.selected[1:3, ]
## coordinates id time Intake.Temperature
## 1 (403540, 5758530) 0 2013-08-11 16:56:06 38
## 21 (403731, 5758060) 20 2013-08-11 16:58:08 32
## 32 (404064, 5758120) 31 2013-08-11 16:59:15 30
## Intake.Pressure CO2 Calculated.MAF Speed Rpm Consumption
## 1 126 0.005875 26.75 4 1339 0.002217
## 21 92 0.005398 24.57 26 1652 0.002037
## 32 112 0.005929 26.99 49 1480 0.002237
spplot(sp.selected, "Consumption", colorkey = TRUE, main = "selected points according to Douglas-Peuker")
Generalization as done above only retains the spatial path, but ignores time and attributes. One could rather easily implement Douglas-Peuker on higher-dimensional curves e.g. \( (x,y,t) \) or \( (x,y,t,Attr_1,Attr_2,...,Attr_n) \) but then one needs a distance metric that combines space, time and attributes – which one would be good?