Using Rayshader for Advanced Geospatial Visualization
R
GIS
Data Science
Author
Pukar Bhandari
Published
January 3, 2023
1 Project Overview
This analysis demonstrates advanced geospatial visualization techniques using R’s rayshader package to create a stunning 3D population density map of Nepal. The project showcases the integration of high-resolution population data with topographic visualization methods, resulting in an interactive and informative representation of demographic patterns across Nepal’s diverse terrain.
# Configure 3D renderingoptions(rgl.useNULL =FALSE)project_crs<-"EPSG:6207"
3 Data Acquisition and Loading
Population Hexagon Data
The analysis utilizes high-resolution population data from Kontur, provided in H3 hexagonal grid format at 400-meter resolution. This dataset offers superior granularity compared to traditional administrative boundary-based population estimates.
# Load population data (H3 hexagons at 400m resolution)nep_hex<-sf::st_read(paste0("/vsigzip/",here::here("posts/nepal-population-density/data/kontur_population_NP_20220630.gpkg.gz")))|>sf::st_transform(project_crs)
Reading layer `population' from data source
`/vsigzip/D:/GitHub/ar-puuk.github.io/posts/nepal-population-density/data/kontur_population_NP_20220630.gpkg.gz'
using driver `GPKG'
Simple feature collection with 122809 features and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 8911286 ymin: 3041698 xmax: 9818046 ymax: 3562066
Projected CRS: WGS 84 / Pseudo-Mercator
Data Sources:
- Population Data: Kontur Population Dataset
- Administrative Boundaries: Kontur Boundaries Dataset
- Coordinate Reference System: EPSG:6207 for accurate representation of Nepal’s geography
4 Data Processing and Visualization Setup
Initial Data Exploration
# Preliminary visualization to verify data qualityggplot2::ggplot(nep_hex)+ggplot2::geom_sf(ggplot2::aes(fill =population), color ="gray66", linewidth =0)+ggplot2::geom_sf(data =nepal_boundary, fill =NA, color ="black", linetype ="dashed", linewidth =1)+ggplot2::theme_minimal()+ggplot2::labs(title ="Nepal Population Distribution", fill ="Population")
Raster Conversion and Matrix Preparation
The conversion from vector hexagon data to raster format is essential for rayshader’s 3D rendering capabilities. This process involves calculating optimal dimensions while preserving the geographic aspect ratio.
# Calculate bounding box and aspect ratiobbox<-sf::st_bbox(nepal_boundary)# Define corner points for aspect ratio calculationbottom_left<-sf::st_point(c(bbox[["xmin"]], bbox[["ymin"]]))|>sf::st_sfc(crs =project_crs)bottom_right<-sf::st_point(c(bbox[["xmax"]], bbox[["ymin"]]))|>sf::st_sfc(crs =project_crs)top_left<-sf::st_point(c(bbox[["xmin"]], bbox[["ymax"]]))|>sf::st_sfc(crs =project_crs)# Calculate dimensionswidth<-sf::st_distance(bottom_left, bottom_right)height<-sf::st_distance(bottom_left, top_left)# Determine aspect ratiosif(width>height){w_ratio<-1h_ratio<-height/width}else{h_ratio<-1w_ratio<-width/height}# Convert to raster with optimal resolutionsize<-1000*2.5pop_raster<-stars::st_rasterize(nep_hex, nx =floor(size*w_ratio), ny =floor(size*h_ratio))# Create matrix for 3D renderingpop_matrix<-matrix(pop_raster$population, nrow =floor(size*w_ratio), ncol =floor(size*h_ratio))
5 Color Scheme Development
MetBrewer Color Palette Selection
The visualization employs the “OKeeffe2” palette from MetBrewer, chosen for its warm earth tones that effectively represent population density while maintaining visual appeal and accessibility.
# Generate color palettecolor<-MetBrewer::met.brewer("OKeeffe2")# Create texture gradient with bias toward higher valuestexture<-grDevices::colorRampPalette(color, bias =3)(256)# Preview color schemecolorspace::swatchplot(texture)
6 3D Visualization Rendering
Interactive 3D Display
# Close any existing RGL windowsrgl::close3d()# Generate 3D visualizationpop_matrix|>rayshader::height_shade(texture =texture)|>rayshader::plot_3d( heightmap =pop_matrix, zscale =250/2.5, solid =FALSE, shadowdepth =0)# Set optimal camera positionrayshader::render_camera( theta =0, phi =30, zoom =0.4, fov =90)
High-Quality Rendering
The final rendering process utilizes rayshader’s advanced lighting system to create publication-ready visualizations with professional-grade quality and resolution.
The final step enhances the visualization with strategic city labels and professional annotation using the {magick} package for image manipulation.
# Load rendered imagepop_raster<-magick::image_read(here::here("posts/nepal-population-density/plots/final_plot.png"))# Define text stylingtext_color<-colorspace::darken(color[7], .5)colorspace::swatchplot(text_color)
# Import and load fontsextrafont::font_import( paths =here::here("posts/nepal-population-density/font"), prompt =FALSE)extrafont::loadfonts(quiet =TRUE)
# Add comprehensive annotationspop_raster|># magick::image_crop(gravity = "center", geometry = "") |> magick::image_annotate("NEPAL", gravity ="northeast", location ="+50+50", color =text_color, size =150,# degrees = 0, font ="Ananda Namaste", weight =700)|>magick::image_annotate("POPULATION DENSITY MAP", gravity ="northeast", location ="+50+200", color =text_color, size =36.5,# degrees = 0, font ="FuturaBT-Medium", weight =500)|># Major urban centersmagick::image_annotate("Kathmandu", gravity ="center", location ="+250-50", color =scales::alpha(text_color, .8), size =30,# degrees = -85, font ="FuturaBT-Medium")|>magick::image_annotate("Pokhara", gravity ="center", location ="-30+35", color =scales::alpha(text_color, .8), size =25,# degrees = -75, font ="FuturaBT-Medium")|>magick::image_annotate("Biratnagar", gravity ="east", location ="+125+100", color =scales::alpha(text_color, .8), size =28,# degrees = -75, font ="FuturaBT-Medium")|># Regional centersmagick::image_annotate("Birgunj", gravity ="center", location ="+130+100", color =scales::alpha(text_color, .8), size =25, font ="FuturaBT-Medium")|>magick::image_annotate("Nepalgunj", gravity ="center", location ="-450+0", color =scales::alpha(text_color, .8), size =24, font ="FuturaBT-Medium")|>magick::image_annotate("Janakpur", gravity ="east", location ="+500+140", color =scales::alpha(text_color, .8), size =22, font ="FuturaBT-Medium")|># Additional urban areasmagick::image_annotate("Itahari", gravity ="east", location ="+230+150", color =scales::alpha(text_color, .8), size =24,# degrees = -75, font ="FuturaBT-Medium")|>magick::image_annotate("Surkhet", gravity ="center", location ="-450-75", color =scales::alpha(text_color, .8), size =20,# degrees = -75, font ="FuturaBT-Medium")|>magick::image_annotate("Tikapur", gravity ="west", location ="+400-125", color =scales::alpha(text_color, .8), size =22,# degrees = -75, font ="FuturaBT-Medium")|>magick::image_annotate("Dhangadhi", gravity ="west", location ="+275-135", color =scales::alpha(text_color, .8), size =18,# degrees = -75, font ="FuturaBT-Medium")|>magick::image_annotate("Bharatpur", gravity ="center", location ="+60+0", color =scales::alpha(text_color, .8), size =20,# degrees = -75, font ="FuturaBT-Medium")|>magick::image_annotate("Butwal", gravity ="center", location ="-125+5", color =scales::alpha(text_color, .8), size =20,# degrees = -75, font ="FuturaBT-Medium")|>magick::image_annotate("Ghorahi\nTulsipur", gravity ="center", location ="-350-100", color =scales::alpha(text_color, .8), size =22,# degrees = -75, font ="FuturaBT-Medium")|># Professional attributionmagick::image_annotate("Visualization by: Pukar Bhandari with Rayshader(@tylermorganwall) | Data: Kontur Population (Released 2022-06-30)", gravity ="southwest", location ="+20+20", color =scales::alpha(text_color, .6),# degrees = 0, font ="FuturaBT-Medium", size =20)|>magick::image_write(here::here("posts/nepal-population-density/nepal_population_density.png"), format ="png", quality =100)
Population Density of Nepal
8 Technical Methodology
Spatial Data Processing Pipeline
Data Transformation: Vector hexagon data converted to standardized coordinate reference system (Nepal Nagarkot TM - ESRI:102306)
Rasterization: Optimal grid resolution calculated to balance detail with computational efficiency
Matrix Conversion: Raster data transformed into height matrix for 3D rendering
Texture Mapping: Population values mapped to color gradients using perceptually uniform color spaces
Rendering Specifications
Resolution: 1980×1080 pixels for high-definition output
Sampling: 300 samples for anti-aliasing and smooth gradients
Lighting: Dual-source lighting system (directional + ambient)
Z-Scale: Optimized vertical exaggeration (250/2.5) for clear population peaks
9 Results and Applications
This methodology produces publication-quality 3D population density visualizations suitable for:
Urban Planning: Identifying population clusters and growth patterns
Infrastructure Development: Targeting high-density areas for transportation networks
Emergency Response: Understanding population distribution for resource allocation
Academic Research: Visual communication of demographic data
10 Conclusion
This workflow demonstrates the power of combining modern R packages for geospatial analysis with advanced 3D visualization techniques. The resulting maps provide intuitive understanding of Nepal’s population distribution patterns, highlighting the concentration of people in the Terai region and major urban centers while showcasing the relatively sparse population in mountainous areas.
The methodology is transferable to other geographic contexts and can be adapted for various demographic or infrastructure planning applications in transportation and urban development projects.