3D bar plot on a map in R

In this post we are going to learn how to plot 3D barplots on a map using R. It is actually not as complicated as it looks. As we will see below, we only need to draw a perspective map using the function persp and then plot the barplots on it.

As usual, we need a shapefile and (for not changing!) we can use this one containing the boundaries of the 38 municipalities of the metropolitan area of Barcelona. We also need a variable containing the information that we want to plot one the map. In our case, we are going to plot random values comprised between 0 and 1. First, we import the shapefile with the readOGR function (package rgdal).

MA=readOGR('MA/MA.shp','MA')
n=dim(MA)[1]

Then we extract the coordinates of the polygons’ centroid (slot labpt). For more information about how to build and manipulate spatial objects with R you can have a look at my previous posts.

cent=cbind(as.numeric(sapply(MA@polygons, function(x) x@labpt[1])),as.numeric(sapply(MA@polygons, function(x) x@labpt[2])))

We generate the random values that we want to plot using the runif function. We obtain one random value per polygon.

myvar=runif(n)

Now, we draw an empty perspective plot using the function persp. The two dimensional vectors x, y and z represents the x-, y- and z-limits. You can adjust the viewing direction with the parameters phi and theta. The function persp returns a projection matrix very useful to superimpose additional graphical elements on the 3D plot and this is exactly what we are going to do.

x=c(minx,maxx)
y=c(miny,maxy)
z=matrix(c(0,0,0,1), 2, 2)

par(mar=c(0,0,0,0))
pmat=persp(x, y, z, xlab = "", ylab = "", zlab="", box=FALSE, border=NA, axes=F, theta=0, phi=20, d=10)

We draw the polygons on our empty perspective plot using the function polygon. The trick here is to use the function trans3D and the projection matrix pmat to project the original coordinates (considered as 3D coordinates) into the 2D plane. It is important to note that in this example each element of the SpatialPolygonDataFrame contains only one Polygon. So we have only one centroid per element, this is not always the case.

for(i in 1:n){
   x=MA[i,1]@polygons[[1]]@Polygons[[1]]@coords[,1]
   y=MA[i,1]@polygons[[1]]@Polygons[[1]]@coords[,2]
   polygon(trans3d(x, y, 0, pmat), col="steelblue3", border="steelblue3") 
} 

Finally, we use the same trick to draw the barplots on the map using the function lines.

for(i in 1:n){
   myz=myvar[i]
   lines(trans3d(rep(cent[i,1], 2), rep(cent[i,2], 2), c(0, myz), pmat), col="#CC6666", lwd=2)
} 

Here is the result:

comments powered by Disqus