Introducing 3D ggplots with Rayshader (R)




Rayshader is a powerful package in R supporting 2D and 3D data visualisation. What amazes me is its details of displaying the 3D plots and it provides lots of customisation. Most importantly, it allows you to directly transform the ggplot2 objects into 3D plot.

Photo by cyda

In this demonstration, the Hong Kong population and housing price data are used. Pre-processed data and the R scripts are uploaded to my Github page. Feel free to download and play around.

Photo by cyda

Package

rayshader

Functionality

rayshader is an open-source package for producing 2D and 3D data visualizations in R. rayshader uses elevation data in a base R matrix and a combination of raytracing, spherical texture mapping, overlays, and ambient occlusion to generate beautiful topographic 2D and 3D maps. In addition to maps, rayshader also allows the user to translate ggplot2 objects into beautiful 3D data visualizations.
The models can be rotated and examined interactively or the camera movement can be scripted to create animations. Scenes can also be rendered using a high-quality pathtracer, rayrender. The user can also create a cinematic depth of field post-processing effect to direct the user’s focus to important regions in the figure. The 3D models can also be exported to a 3D-printable format with a built-in STL export function.

Demonstration

  1. Generate a 2D map plot with ggplot2
  2. Transform the 2D plot into a 3D plot with rayshader

Data

Population Data.xlsx

Photo by cyda

real_estate_master_df.csv

Photo by cyda


Task 1: Generate a 2D map plot with ggplot2
library(sp)
hkmap = readRDS("HKG_adm1.rds") # geo data of HK map
# Preprocessing
map_data = data.frame(id=hkmap$ID_1, Code=hkmap$HASC_1, Eng_name=hkmap$NAME_1)
map_data$Code = gsub('HK.', '', as.character(map_data$Code))
map_data = merge(map_data, district_name, by = 'Eng_name')
hkmapdf = fortify(hkmap)
map_data = merge(hkmapdf, map_data, by="id")
map_data = merge(map_data, population, by = "Chi_name")
map_data$Population = as.numeric(map_data$Population)
library(ggplot2)
# Map
map_bg = ggplot(map_data, aes(long, lat, group=group, fill = Population)) +
  geom_polygon() + # Shape
  scale_fill_gradient(limits=range(map_data$Population), 
                      low="#FFF3B0", high="#E09F3E") + 
  layer(geom="path", stat="identity", position="identity", 
       mapping=aes(x=long, y=lat, group=group, 
                   color=I('#FFFFFF'))) + 
  theme(legend.position = "none", 
        axis.line=element_blank(), 
        axis.text.x=element_blank(), axis.title.x=element_blank(),
        axis.text.y=element_blank(), axis.title.y=element_blank(),
        axis.ticks=element_blank(), 
        panel.background = element_blank()) # Clean Everything
map_bg
# Save as PNG
xlim = ggplot_build(map_bg)$layout$panel_scales_x[[1]]$range$range
ylim = ggplot_build(map_bg)$layout$panel_scales_y[[1]]$range$range
ggsave('map_bg1.png', width = diff(xlim)*100, height = diff(ylim)*100, units = "cm")
Step-by-step explanation:
hkmap = readRDS("HKG_adm1.rds") # geo data of HK map
To plot a map of any regions, you need to have the location data on hand. Please download the rds file from the Github repo first.
library(ggplot2)
# Map
map_bg = ggplot(map_data, aes(long, lat, group=group, fill = Population)) +
  geom_polygon() + # Shape
  scale_fill_gradient(limits=range(map_data$Population), 
                      low="#FFF3B0", high="#E09F3E") + 
  layer(geom="path", stat="identity", position="identity", 
       mapping=aes(x=long, y=lat, group=group, 
                   color=I('#FFFFFF')))
ggplot2 is another wonderful package in R providing many kinds of charts like bar chart, pie chart, heat map … you name it. In the script shown above, there are 2 main things done. On the line, scale_fill_gradient, the polygons are shown in the map are filled by colour based on population data. Regions with darker colour mean that they have a higher population. Then, layer gives a borderline to each region which is the white line shown.

Photo by cyda

map_bg + 
  theme(legend.position = "none", 
        axis.line=element_blank(), 
        axis.text.x=element_blank(), axis.title.x=element_blank(),
        axis.text.y=element_blank(), axis.title.y=element_blank(),
        axis.ticks=element_blank(), 
        panel.background = element_blank()) # Clean Everything
map_bg
The reason for creating this plot is to prepare a background image for the 3D plot that we are going to work out in task 2. Thus, I remove all the gridlines, x-axis and y-axis by using theme.

Output: map_bg.png

# Save as PNG
xlim = ggplot_build(map_bg)$layout$panel_scales_x[[1]]$range$range
ylim = ggplot_build(map_bg)$layout$panel_scales_y[[1]]$range$range
ggsave(‘map_bg.png’, width = diff(xlim)*40, height = diff(ylim)*40, units = “cm”)
Lastly, ggsave is used to export the map_bg with the right scale derived from the plot coordinates.

Task 2: Transform the 2D plot into a 3D plot with rayshader
Data of real_estate_master_df.csv are scraped from Hong Kong Property.
# 2D Plot
library(ggplot2)
library(grid)
estate_price = ggplot(estate_df) + 
  annotation_custom(rasterGrob(hk_map_bg, width=unit(1,"npc"), 
                               height=unit(1,"npc")), 
                    -Inf, Inf, -Inf, Inf) + 
  xlim(xlim[1],xlim[2]) + # x-axis Mapping
  ylim(ylim[1],ylim[2]) + # y-axis Mapping
  geom_point(aes(x=Longitude, y=Latitude, color=apr_price),size=2) + 
  scale_colour_gradient(name = '成交呎價(實)\n(HKD)', 
                        limits=range(estate_df$apr_price), 
                        low="#FCB9B2", high="#B23A48") + 
  theme(axis.line=element_blank(), 
        axis.text.x=element_blank(), axis.title.x=element_blank(),
        axis.text.y=element_blank(), axis.title.y=element_blank(),
        axis.ticks=element_blank(), 
        panel.background = element_blank()) # Clean Everything
estate_price
Step-by-step explanation:
  annotation_custom(rasterGrob(hk_map_bg, width=unit(1,"npc"), 
                               height=unit(1,"npc")), 
                    -Inf, Inf, -Inf, Inf) + 
  xlim(xlim[1],xlim[2]) + # x-axis Mapping
  ylim(ylim[1],ylim[2]) + # y-axis Mapping
annotation_custom is used to import the prepared map background image to the plot. xlim(xlim[1],xlim[2]) and ylim(ylim[1],ylim[2]) fix the x- and y-axis of the plot and make them the same as the map_bg.
geom_point(aes(x=Longitude, y=Latitude, color=apr_price),size=2) + 
  scale_colour_gradient(name = '成交呎價(實)\n(HKD)', 
                        limits=range(estate_df$apr_price), 
                        low="#FCB9B2", high="#B23A48") +
geom_point adds the data points of estate price to the plot.
scale_colour_gradient fills the data points with color based on the variable, “apr_price”, please note that the function used here is different from the one used in task 1. (scale_colour_gradient for geom_point, scale_fill_gradient for geom_polygon)

Photo by cyda

# 3D Plot
library(rayshader)
plot_gg(estate_price, multicore = TRUE, width = diff(xlim)*10 ,height=diff(ylim)*10, fov = 70, scale = 300)
# Close Windows
rgl.close()
The script to convert the 2D plot to a 3D plot is very simple. The 3D plot is shown in the rgl window. However, remember to close the window after viewing the plot. As the output of the rgl window will not be clear, so the next plot will be overlapped with the previous plot.

Photo by cyda

Hope you enjoy my article =)
To know more about using rayshader, please go to the official website of rayshader.

Comments

Post a Comment