Z - Values for polygon (shapefile) in R

Question

my goal is to create a 3D-Visualization in R. I have a shapefile of urban districts (Ortsteile) in Berlin and want to highlight the value (inhabitants/km²) as a z-value. I have implemented the shapefile into R and coloured the value for desnity ("Einwohnerd") as followed:

library(rgdal)
library(sp)

berlin=readOGR(dsn="C...etc.", layer="Ortsteile")

berlin@data

col <- rainbow(length(levels(berlin@data$Name)))
spplot(berlin, "Einwohnerd", col.regions=col, main="Ortsteil Berlins", sub="Datensatz der Stadt Berlin", lwd=.8, col="black")

How it is posible to refer a certain polygon (urban district) to a z-value (inhabitant/km²) and how can I highlight this z-value?

Hope that someone will have an answer! Best regars SB

Thanks for the answer, but I am still on my wy to find out the best to use the density as z-value so that I can create a 3D Model. I found out that it is not possible to use the polygons of the shape but that it is possible to rasterize the polygon and to use a matrix for a different perspective and rotation.

Here is the code but the final 3D visualization looks not sharp and good enough. Maybe it would be better to calculate the the z-value in anther way so that the first values did not start so high or to use the center of the polygon and than to draw a column in z-direction:

library(rgdal)
library(sp)

setwd("C:\\...")
berlin=readOGR(dsn="C:\\...\\Ortsteile", layer="Ortsteile") 

col <- rainbow(length(levels(berlin@data$Name)))  
spplot(berlin, "Einwohnerd", col.regions=col, main="Ortsteil Berlins",                 
sub="Datensatz    der Stadt Berlin", lwd=.8, col="black")

library(raster)

raster <- raster(nrows=100, ncols=200, extent(berlin)) 

test <- rasterize(berlin, raster, field="Einwohnerd")

persp(test, theta = 40, phi = 40, col = "gold", border = NA, shade = 0.5)  

for(i in seq(0,90,10)){     
persp(test, theta = 40, phi = i, col = "gold", border = NA, shade = 0.5)
}

library(rgl)         
library(colorRamps)
mat <- matrix(test[], nrow=test@nrows, byrow=TRUE)
image(mat)
persp3d(z = mat, clab = "m")
persp3d(z = mat, col = rainbow(10),border = "black")
persp3d(z = mat, facets = FALSE, curtain = TRUE)

Answer 1

Is this what you had in mind?

library(ggplot2)
library(rgdal)           # for readOGR(...) and spTransform(...)
library(RColorBrewer)    # for brewer.pal(...)

setwd("<directory with shapefile>")
map <- readOGR(dsn=".",layer="Ortsteile")
map <- spTransform(map,CRS=CRS("+init=epsg:4839"))
map.data <- data.frame(id=rownames(map@data), map@data)
map.df   <- fortify(map)
map.df   <- merge(map.df,map.data,by="id")
ggplot(map.df, aes(x=long, y=lat, group=group))+
  geom_polygon(aes(fill=Einwohnerd))+
  geom_path(colour="grey")+
  scale_fill_gradientn(colours=rev(brewer.pal(10,"Spectral")))+
  theme(axis.text=element_blank())+
  labs(title="Berlin Ortsteile", x="", y="")+
  coord_fixed()

Explanation

This is a great question, in that it provides an example of a very basic choropleth map using ggplot in R.

Shapefiles can be read into R using readOGR(...), producing SpatialDataFrame objects. The latter have basically two sections: a polygons section containing the coordinates of the polygon boundaries, and a data section containing information from the attributes table in the shapefile. These can be referenced, respectively, as map@polygons and map@data.

The code above reads the shapefile and transforms the coordinates to epsg:4839. Then we prepend the polygon ids (stored in the rownames) to the other information in map@data, creating map.data. Then we use the fortify(...) function in ggplot to convert the polygons to a dataframe suitable for plotting (map.df). This dataframe has a column id which corresponds to the id column in map.data. Then we merge the attribute information (map.data) into map.df based on the id column.

The ggplot calls create the map layers and render the map, as follows:

ggplot:       set the default dataset to map.df; identify x- and y-axis columns
geom_polygon: identify column for fill (color of polygon)
geom_path:    polygon boundaries
theme:        turn off axis text
labs:         title, turn off x- and y-axis labels
coord_fixed:  ensures that the map is not distorted

A note on scale_fill_gradientn(...): this function assigns colors to the fill values by interpolating a color palette provided in the colours= parameter. Here we use the Spectral palette from www.colorbrewer.org. Unfotrunately, this palette has the colors revered (blue - red), so we use rev(...) to reverse the color order (high=red, low=blue). If you prefer the more highly saturated colors common in matlab, use library(colorRamps) and replace the call to scale_fill_gradientn(...) with:

  scale_fill_gradientn(colours=matlab.like(10))+