**everyday analytics**, and kindly contributed to R-bloggers)

A little while back there was an article in blogTO about how a reddit user had used data from Toronto’s Open Data initiative to produce a rather cool-looking map of all the locations of all the traffic signals here in the city.

It’s neat because as the author on blogTO notes, it is recognizable as Toronto without any other geographic data being plotted – the structure of the city comes out in the data alone.

Still, I thought it’d be interesting to see as a geographic heat map, and also a good excuse to fool around with mapping using Rgooglemaps.

The finished product below:

Despite my best efforts with transparency (using my helper function), it’s difficult for anything but the city core to really come out in the intensity map.

The image without the Google maps tile, and the coordinates rotated, shows the density a little better in the green-yellow areas:

And it’s also straightforward to produce a duplication of the original black and white figure:

The R code is below. Interpolation is using the trusty kde2d function from the MASS library and a rotation is applied for the latter two figures, so that the grid of Toronto’s streets faces ‘up’ as in the original map.

# Toronto Traffic Signals Heat Map

# Myles Harrison

# http://www.everydayanalytics.ca

# Data from Toronto Open Data Portal:

# http://www.toronto.ca/open

library(MASS)

library(RgoogleMaps)

library(RColorBrewer)

source('colorRampPaletteAlpha.R')

# Read in the data

data <- read.csv(file="traffic_signals.csv", skip=1, header=T, stringsAsFactors=F)

# Keep the lon and lat data

rawdata <- data.frame(as.numeric(data$Longitude), as.numeric(data$Latitude))

names(rawdata) <- c("lon", "lat")

data <- as.matrix(rawdata)

# Rotate the lat-lon coordinates using a rotation matrix

# Trial and error lead to pi/15.0 = 12 degrees

theta = pi/15.0

m = matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), nrow=2)

data <- as.matrix(data) %*% m

# Reproduce William's original map

par(bg='black')

plot(data, cex=0.1, col="white", pch=16)

# Create heatmap with kde2d and overplot

k <- kde2d(data[,1], data[,2], n=500)

# Intensity from green to red

cols <- rev(colorRampPalette(brewer.pal(8, 'RdYlGn'))(100))

par(bg='white')

image(k, col=cols, xaxt='n', yaxt='n')

points(data, cex=0.1, pch=16)

# Mapping via RgoogleMaps

# Find map center and get map

center <- rev(sapply(rawdata, mean))

map <- GetMap(center=center, zoom=11)

# Translate original data

coords <- LatLon2XY.centered(map, rawdata$lat, rawdata$lon, 11)

coords <- data.frame(coords)

# Rerun heatmap

k2 <- kde2d(coords$newX, coords$newY, n=500)

# Create exponential transparency vector and add

alpha <- seq.int(0.5, 0.95, length.out=100)

alpha <- exp(alpha^6-1)

cols2 <- addalpha(cols, alpha)

# Plot

PlotOnStaticMap(map)

image(k2, col=cols2, add=T)

points(coords$newX, coords$newY, pch=16, cex=0.3)

This a neat little start and you can see how this type of thing could easily be extended to create a generalized mapping tool, stood up as a web service for example (they’re out there). Case in point: Google Fusion Tables. I’m unsure as to what algorithm they use but I find it less satisfying, looks like some kind of simple point blending:

As always, all the code is on github.

**leave a comment**for the author, please follow the link and comment on his blog:

**everyday analytics**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...