Using {Rayshaders} package to visualize 3D map in R

What is the {Rayshaders} package?

“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, hillshading algorithms, and overlays to generate stunning 2D and 3D maps. In addition to maps, rayshader also allows the user to translate ggplot2 objects into beautiful 3D data visualizations.” https://www.rayshader.com/

Here are the instructions for creating a simple 3D map for Connecticut’s median household income by census tract.

Set working directory.

setwd("~/Documents/Visualization 3D map Using Rayshaders Package in R")

View files in working directory.

list.files()

##  [1] "2022_Disproportionately_Impacted_Areas_20240212.csv"     
##  [2] "tl_2019_09_tract.cpg"                                    
##  [3] "tl_2019_09_tract.dbf"                                    
##  [4] "tl_2019_09_tract.prj"                                    
##  [5] "tl_2019_09_tract.shp"                                    
##  [6] "tl_2019_09_tract.shp.ea.iso.xml"                         
##  [7] "tl_2019_09_tract.shp.iso.xml"                            
##  [8] "tl_2019_09_tract.shx"                                    
##  [9] "Using {Rayshaders} package to visualize 3D map in R.Rmd" 
## [10] "Using-{Rayshaders}-Package-to-visualize-3D-map-in-R.html"
## [11] "Using-{Rayshaders}-package-to-visualize-3D-map-in-R.Rmd" 
## [12] "Visualization 3D map Using Rayshaders Package in R.R"

Upload libraries.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(ggplot2)
library(dplyr)
library(rayshader)
library(viridis)

## Loading required package: viridisLite

library(viridisLite)
library(magick)

## Linking to ImageMagick 6.9.12.93
## Enabled features: cairo, fontconfig, freetype, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fftw, ghostscript, x11

library(av)
library(janitor)

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(plotly)

## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

Upload data and map files.

# Import data file
mydata<- read_csv("2022_Disproportionately_Impacted_Areas_20240212.csv")

## Rows: 883 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): GEOID, Town(s)
## dbl (5): Median Household Income, Population, Conviction Count, Conviction R...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

print(mydata)

## # A tibble: 883 × 7
##    GEOID       `Town(s)`    Median Household Inc…¹ Population `Conviction Count`
##    <chr>       <chr>                         <dbl>      <dbl>              <dbl>
##  1 09003980001 East Granby                      NA          0                0  
##  2 09003510100 East Hartfo…                  51089       1873              154. 
##  3 09003510900 East Hartfo…                  96630       3684              114. 
##  4 09003511400 East Hartfo…                  84918       2167               97  
##  5 09003514101 Manchester                    84853       3096               66.1
##  6 09003514800 Manchester                    50511       3114              480. 
##  7 09001044600 Norwalk                      234152       3583               32.1
##  8 09001045101 Wilton                       209899       4018               21.6
##  9 09001045102 Wilton                       189479       6027               62.6
## 10 09003511300 East Hartfo…                  42527       3503              389  
## # ℹ 873 more rows
## # ℹ abbreviated name: ¹`Median Household Income`
## # ℹ 2 more variables: `Conviction Rate` <dbl>,
## #   `Disproportionately Impacted Area (DIA)` <dbl>

# Import shapefile
mymap<- st_read("tl_2019_09_tract.shp")

## Reading layer `tl_2019_09_tract' from data source 
##   `/Users/ramonrodriguez-santana/Documents/Visualization 3D map Using Rayshaders Package in R/tl_2019_09_tract.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 833 features and 12 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.72777 ymin: 40.95094 xmax: -71.78724 ymax: 42.05051
## Geodetic CRS:  NAD83

print(mymap)

## Simple feature collection with 833 features and 12 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.72777 ymin: 40.95094 xmax: -71.78724 ymax: 42.05051
## Geodetic CRS:  NAD83
## First 10 features:
##    STATEFP COUNTYFP TRACTCE       GEOID NAME          NAMELSAD MTFCC FUNCSTAT
## 1       09      003  405700 09003405700 4057 Census Tract 4057 G5020        S
## 2       09      003  510200 09003510200 5102 Census Tract 5102 G5020        S
## 3       09      003  510400 09003510400 5104 Census Tract 5104 G5020        S
## 4       09      015  901100 09015901100 9011 Census Tract 9011 G5020        S
## 5       09      003  510300 09003510300 5103 Census Tract 5103 G5020        S
## 6       09      003  492400 09003492400 4924 Census Tract 4924 G5020        S
## 7       09      003  405200 09003405200 4052 Census Tract 4052 G5020        S
## 8       09      003  415300 09003415300 4153 Census Tract 4153 G5020        S
## 9       09      003  492500 09003492500 4925 Census Tract 4925 G5020        S
## 10      09      003  494100 09003494100 4941 Census Tract 4941 G5020        S
##        ALAND  AWATER    INTPTLAT     INTPTLON                       geometry
## 1    3590069   77261 +41.6767836 -072.9749700 POLYGON ((-72.99138 41.6799...
## 2    2897973  309925 +41.7666027 -072.6540100 POLYGON ((-72.66507 41.7665...
## 3    2709322   87828 +41.7716475 -072.6299729 POLYGON ((-72.64662 41.7642...
## 4  157165208 3070895 +41.9662084 -072.0227864 POLYGON ((-72.10217 42.0288...
## 5    2388034       0 +41.7858322 -072.6246204 POLYGON ((-72.6389 41.77622...
## 6    3008123    3059 +41.7031581 -072.6896542 POLYGON ((-72.70626 41.7029...
## 7    3629846       0 +41.6796950 -072.9124225 POLYGON ((-72.9269 41.68972...
## 8     541180       0 +41.6673973 -072.7681828 POLYGON ((-72.7723 41.66874...
## 9    3223856       0 +41.6954213 -072.6880119 POLYGON ((-72.70527 41.6980...
## 10   5044765       0 +41.6766129 -072.7164631 POLYGON ((-72.73047 41.6838...

Join mymap and mydata files by GEOID, using the inner_join() function.

# Join map and data files by GEOID
map_and_data_shp<- inner_join(mymap, mydata, by ="GEOID")
print(map_and_data_shp)

## Simple feature collection with 773 features and 18 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.72777 ymin: 40.95094 xmax: -71.78724 ymax: 42.05051
## Geodetic CRS:  NAD83
## First 10 features:
##    STATEFP COUNTYFP TRACTCE       GEOID    NAME             NAMELSAD MTFCC
## 1       09      003  405700 09003405700    4057    Census Tract 4057 G5020
## 2       09      003  510200 09003510200    5102    Census Tract 5102 G5020
## 3       09      003  510400 09003510400    5104    Census Tract 5104 G5020
## 4       09      003  510300 09003510300    5103    Census Tract 5103 G5020
## 5       09      003  492400 09003492400    4924    Census Tract 4924 G5020
## 6       09      003  405200 09003405200    4052    Census Tract 4052 G5020
## 7       09      003  415300 09003415300    4153    Census Tract 4153 G5020
## 8       09      003  492500 09003492500    4925    Census Tract 4925 G5020
## 9       09      003  494100 09003494100    4941    Census Tract 4941 G5020
## 10      09      003  494201 09003494201 4942.01 Census Tract 4942.01 G5020
##    FUNCSTAT   ALAND AWATER    INTPTLAT     INTPTLON       Town(s)
## 1         S 3590069  77261 +41.6767836 -072.9749700       Bristol
## 2         S 2897973 309925 +41.7666027 -072.6540100 East Hartford
## 3         S 2709322  87828 +41.7716475 -072.6299729 East Hartford
## 4         S 2388034      0 +41.7858322 -072.6246204 East Hartford
## 5         S 3008123   3059 +41.7031581 -072.6896542  Wethersfield
## 6         S 3629846      0 +41.6796950 -072.9124225       Bristol
## 7         S  541180      0 +41.6673973 -072.7681828   New Britain
## 8         S 3223856      0 +41.6954213 -072.6880119  Wethersfield
## 9         S 5044765      0 +41.6766129 -072.7164631     Newington
## 10        S 6059664      0 +41.6611974 -072.7392738     Newington
##    Median Household Income Population Conviction Count Conviction Rate
## 1                    38650       2075            292.1            0.14
## 2                    56103       2578            364.2            0.14
## 3                    44929       6535            999.8            0.15
## 4                    49242       4087            567.5            0.14
## 5                    90735       2866            117.3            0.04
## 6                    65156       4201            239.2            0.06
## 7                    34963       2067            578.9            0.28
## 8                   113871       3345             35.6            0.01
## 9                    82750       5965            208.0            0.03
## 10                   77536       4617            137.3            0.03
##    Disproportionately Impacted Area (DIA)                       geometry
## 1                                       1 POLYGON ((-72.99138 41.6799...
## 2                                       1 POLYGON ((-72.66507 41.7665...
## 3                                       1 POLYGON ((-72.64662 41.7642...
## 4                                       1 POLYGON ((-72.6389 41.77622...
## 5                                       0 POLYGON ((-72.70626 41.7029...
## 6                                       0 POLYGON ((-72.9269 41.68972...
## 7                                       1 POLYGON ((-72.7723 41.66874...
## 8                                       0 POLYGON ((-72.70527 41.6980...
## 9                                       0 POLYGON ((-72.73047 41.6838...
## 10                                      0 POLYGON ((-72.75431 41.6721...

Clean map_and_data_shp file variables using the clean_names() function and view new renamed variables.

map_and_data_shp <- map_and_data_shp %>%
  clean_names()
# List renamed variable
ls(map_and_data_shp)

##  [1] "aland"                               
##  [2] "awater"                              
##  [3] "conviction_count"                    
##  [4] "conviction_rate"                     
##  [5] "countyfp"                            
##  [6] "disproportionately_impacted_area_dia"
##  [7] "funcstat"                            
##  [8] "geoid"                               
##  [9] "geometry"                            
## [10] "intptlat"                            
## [11] "intptlon"                            
## [12] "median_household_income"             
## [13] "mtfcc"                               
## [14] "name"                                
## [15] "namelsad"                            
## [16] "population"                          
## [17] "statefp"                             
## [18] "town_s"                              
## [19] "tractce"

View interactive map using ggplotly() function.

p<- ggplot(map_and_data_shp) +
  geom_sf(aes(fill = median_household_income)) +
  scale_fill_gradient(low = "#56B1F7", high = "#132B43") +
  ggtitle("Connecticut's MHI by Census Tract (2022)") 

ggplotly(p)

Transform coordinate reference system (crs) to WGS 84 format by using crs = 4326.

ct_mhi_shp<- st_transform(map_and_data_shp, crs = 4326)
print(ct_mhi_shp)

## Simple feature collection with 773 features and 18 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.72777 ymin: 40.95094 xmax: -71.78724 ymax: 42.05051
## Geodetic CRS:  WGS 84
## First 10 features:
##    statefp countyfp tractce       geoid    name             namelsad mtfcc
## 1       09      003  405700 09003405700    4057    Census Tract 4057 G5020
## 2       09      003  510200 09003510200    5102    Census Tract 5102 G5020
## 3       09      003  510400 09003510400    5104    Census Tract 5104 G5020
## 4       09      003  510300 09003510300    5103    Census Tract 5103 G5020
## 5       09      003  492400 09003492400    4924    Census Tract 4924 G5020
## 6       09      003  405200 09003405200    4052    Census Tract 4052 G5020
## 7       09      003  415300 09003415300    4153    Census Tract 4153 G5020
## 8       09      003  492500 09003492500    4925    Census Tract 4925 G5020
## 9       09      003  494100 09003494100    4941    Census Tract 4941 G5020
## 10      09      003  494201 09003494201 4942.01 Census Tract 4942.01 G5020
##    funcstat   aland awater    intptlat     intptlon        town_s
## 1         S 3590069  77261 +41.6767836 -072.9749700       Bristol
## 2         S 2897973 309925 +41.7666027 -072.6540100 East Hartford
## 3         S 2709322  87828 +41.7716475 -072.6299729 East Hartford
## 4         S 2388034      0 +41.7858322 -072.6246204 East Hartford
## 5         S 3008123   3059 +41.7031581 -072.6896542  Wethersfield
## 6         S 3629846      0 +41.6796950 -072.9124225       Bristol
## 7         S  541180      0 +41.6673973 -072.7681828   New Britain
## 8         S 3223856      0 +41.6954213 -072.6880119  Wethersfield
## 9         S 5044765      0 +41.6766129 -072.7164631     Newington
## 10        S 6059664      0 +41.6611974 -072.7392738     Newington
##    median_household_income population conviction_count conviction_rate
## 1                    38650       2075            292.1            0.14
## 2                    56103       2578            364.2            0.14
## 3                    44929       6535            999.8            0.15
## 4                    49242       4087            567.5            0.14
## 5                    90735       2866            117.3            0.04
## 6                    65156       4201            239.2            0.06
## 7                    34963       2067            578.9            0.28
## 8                   113871       3345             35.6            0.01
## 9                    82750       5965            208.0            0.03
## 10                   77536       4617            137.3            0.03
##    disproportionately_impacted_area_dia                       geometry
## 1                                     1 POLYGON ((-72.99138 41.6799...
## 2                                     1 POLYGON ((-72.66507 41.7665...
## 3                                     1 POLYGON ((-72.64662 41.7642...
## 4                                     1 POLYGON ((-72.6389 41.77622...
## 5                                     0 POLYGON ((-72.70626 41.7029...
## 6                                     0 POLYGON ((-72.9269 41.68972...
## 7                                     1 POLYGON ((-72.7723 41.66874...
## 8                                     0 POLYGON ((-72.70527 41.6980...
## 9                                     0 POLYGON ((-72.73047 41.6838...
## 10                                    0 POLYGON ((-72.75431 41.6721...

Rename median_household_income to ‘ct_mhi’ and create new ct_mhi_shp shapefile .

ct_mhi_shp = mutate(ct_mhi_shp, ct_mhi = median_household_income)

Visualization of 3D map using {rayshader} package

Create 3D Map visualization.

ggCTmhi <- ggplot(data=ct_mhi_shp) +
  geom_sf(aes(fill = ct_mhi)) +
  scale_fill_viridis() +
  ggtitle("Connecticut's MHI by Census Tract (2022)") +
  theme_bw()

Plot 3D Map visualization.
Note. Although you cannot view the 3D map in this presentation, you will be able to view it in your computer.

plot_gg(ggCTmhi,multicore = TRUE, width = 5, height = 5, scale = 200,
        windowsize=c(1280,720), zoom = 0.60, phi = 50, sunangle = 120, theta = 45)

Finally, spacial thanks to the authors/developers of the {Rayshaders} package.

For more information regarding the package please visit:
https://www.rayshader.com/

A.M.D.G.

Using {Rayshaders} package to visualize 3D map in R

Ramon Rodriguez-Santana, MBA, MPH

2024-02-13

What is the {Rayshaders} package?

Visualization of 3D map using {rayshader} package