public
Last active

Plotting 3D maps using OpenStreetMap and RGL.

  • Download Gist
3d_mapping
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
# Plotting 3D maps using OpenStreetMap and RGL. For info see:
# http://geotheory.co.uk/blog/2013/04/26/plotting-3d-maps-with-rgl/
 
map3d <- function(map, ...){
if(length(map$tiles)!=1){stop("multiple tiles not implemented") }
nx = map$tiles[[1]]$xres
ny = map$tiles[[1]]$yres
xmin = map$tiles[[1]]$bbox$p1[1]
xmax = map$tiles[[1]]$bbox$p2[1]
ymin = map$tiles[[1]]$bbox$p1[2]
ymax = map$tiles[[1]]$bbox$p2[2]
xc = seq(xmin,xmax,len=ny)
yc = seq(ymin,ymax,len=nx)
colours = matrix(map$tiles[[1]]$colorData,ny,nx)
m = matrix(0,ny,nx)
surface3d(xc,yc,m,col=colours, ...)
}
 
#Sys.setenv(NOAWT=1) # Mac users: fixes an OSM/X11 issue that may arise
require(rgl)
require(OpenStreetMap)
require(ggplot2)
require(maptools)
 
# download map tile (the '8' parameter for map resolution)
lat <- c(51.7, 51.3); lon <- c(-0.53, 0.3)
map <- openproj(openmap(c(lat[1],lon[1]),c(lat[2],lon[2]), 8, 'osm'))
 
# import London rents data (originally from London Data Store)
rents <- read.csv("http://bit.ly/YZRYEC", header=T)
head(rents)
 
# create xyz matrix for line heights (row pairs for segment points)
m <- mat.or.vec(66, 3)
for(i in 1:66) for(j in 1:3) m[i,j] = rents[ceiling(i/2),j+2]
for(i in 1:33) m[i*2,3] = 0; m[,3] = m[,3]/15000
head(m)
 
# draw map
run <- function(){
open3d(windowRect=c(100,100,800,600))
map3d(map, lit=F)
segments3d(m, lwd=5, col='red', alpha=0.6)
play3d(spin3d(axis=c(0,0,0.5), rpm=24), duration=2.5)
}
run()
 
# adjust to custom view and save to file
rgl.snapshot("output.png", fmt="png", top=TRUE)

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.