Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Plotting 3D maps using OpenStreetMap and RGL.
# Plotting 3D maps using OpenStreetMap and RGL. For info see:
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
# 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("", header=T)
# 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
# draw map
run <- function(){
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)
# adjust to custom view and save to file
rgl.snapshot("output.png", fmt="png", top=TRUE)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment