<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/create-3d-county-maps-using-density-as-z-axis/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://rlang.io/create-3d-county-maps-using-density-as-z-axis/"> RLang.io | R Language Programming</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This is going to be a bit longer than some of my previous tutorials as it covers a walkthrough for sourcing data, scraping tables, cleaning, and generating the 3D view below which you can springboard from with the help of the <code>rgl</code> package. The heavy lifting is done with <code>ggplot</code> and <code>rayshader</code>.</p>
<p><img alt="Rendered Visualization of MN Election Results in 2018 for Governor" title="Rendered Visualization of MN Election Results in 2018 for Governor" src="https://i2.wp.com/rlang.io/wp-content/uploads/2018/11/MN-Election-Results-rgl.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h1>Rayshader</h1>
<blockquote>
<p>
rayshader is an open source R package for producing 2D and 3D hillshaded maps of elevation matrices using a combination of raytracing, spherical texture mapping, and ambient occlusion.
</p>
</blockquote>
<p><a rel="nofollow" href="https://github.com/tylermorganwall/rayshader" target="_blank">Rayshader on GitHub</a></p>
<p>This amazing package is what inspired this tutorial. I saw one of the author’s tweets, remembered the <a rel="nofollow" href="https://blueshift.io/election-2016-county-map.html" target="_blank">3D map from BlueShift</a> I saw from the 2016 election results.</p>
<p>One thing to keep in mind is this will not be as fast as some other solutions (plotly for example), and the interactive element will <em>likely</em> not be shareable. There is a function called <code>writeWebGL</code> within rgl that I have not had much luck with, but it does exist. Rasyshader allows for creating programatic flyovers and animations though.</p>
<h1>The R Script</h1>
<h2>Including packages</h2>
<p>These are the packages I used for the majority of the code. I will include the ones for creating rasters and rayshading later in this write-up.</p>

<pre>library(ggplot2)
library(ggmap)
library(maps)
library(mapdata)
library(stringr)
library(dplyr)
library(rvest)
library(magrittr)
library(ggthemes)</pre>

<h2>County Information</h2>
<p>This tutorial is for Minnesota, but other regions can be used without changing too much. The biggest difference I noticed was abbreviations for county names.</p>

<pre>usa &lt;- map_data(&quot;usa&quot;)
states &lt;- map_data(&quot;state&quot;)
mn_df &lt;- subset(states, region == &quot;minnesota&quot;)
counties &lt;- map_data(&quot;county&quot;)
mn_county &lt;- subset(counties, region == &quot;minnesota&quot;)
mn_county$pos &lt;- 1:nrow(mn_county)</pre>

<p>The dataframe <code>mn_county</code> will be reused as the primary dataframe in this tutorial, and will have information from difference sources appended to this dataframe.</p>
<h2>Population Information</h2>

<pre>webpage &lt;- read_html(&quot;https://en.wikipedia.org/wiki/List_of_counties_in_Minnesota&quot;)
tbls &lt;- html_nodes(webpage, &quot;table&quot;)
# Shows all tables
tbls</pre>

<p>We only need one table. To determine which table we print out the table. It is the one with sortable in the class name, but there are other ways of determining which one. The results from <code>tbls</code> will give something that looks like this…</p>

<pre>{xml_nodeset (3)}
[1] <table>

<p>To make this easier I just chose to use the first table to read, and used the first element (dataframe) from the <code>html_table</code> return.</p>

<pre>wiki &lt;- html_table(tbls[1],fill = T)[[1]]
# Remove Citations in Column Names
names(wiki) &lt;- gsub(&quot;\\[.*$&quot;,&quot;&quot;,names(wiki))
# Convert to Numeric and Remove Weird Characters
wiki$Population &lt;- wiki$Population %>% 
  gsub(&quot;^.*<img alt="♠" src="https://i2.wp.com/s.w.org/images/core/emoji/11/72x72/2660.png?w=456&ssl=1" class="wp-smiley" data-recalc-dims="1">&quot;,&quot;&quot;,.) %>%
  gsub(&quot;[^0-9\\.]&quot;,&quot;&quot;,.) %>%
  as.numeric
# Convert to Numeric and Remove Weird Characters
wiki$Area &lt;- wiki$Area %>% 
  gsub(&quot;^[0-9]+<img alt="♠" src="https://i2.wp.com/s.w.org/images/core/emoji/11/72x72/2660.png?w=456&ssl=1" class="wp-smiley" data-recalc-dims="1">&quot;,&quot;&quot;,.) %>%
  gsub(&quot;sq.*$&quot;,&quot;&quot;,.) %>%
  gsub(&quot;[^0-9\\.]&quot;,&quot;&quot;,.) %>%
  as.numeric
# Column not needed
wiki$Map &lt;- NULL
# Remove &quot; County&quot; from County Names
# One off replacement for &quot;saint louis&quot;
wiki$County &lt;- gsub(&quot;( County|\\.)&quot;,&quot;&quot;,wiki$County) %>% 
  tolower %>% gsub(&quot;saint louis&quot;,&quot;st louis&quot;,.)
# Just makes it easier to merge later
names(wiki)[1] &lt;- &quot;subregion&quot;

# Append to mn_county
mn_county &lt;- merge(mn_county,wiki,by=&quot;subregion&quot;,all.x=T)
mn_county$density &lt;- mn_county$Population / mn_county$Area
# Unused
mn_county$bin &lt;- NULL</pre>

<p>The county density is used for our <em>“elevation”</em>.</p>
<p><img alt="MN Density" title="MN Density" src="https://i1.wp.com/rlang.io/wp-content/uploads/2018/11/MN-Density.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Obviously this isn’t really a worthwhile result yet, but it is getting there. The code for the image above is essentially the same as the code near the end, but with the color layer ommited.</p>
<h2>Governor Information</h2>

<pre>webpage &lt;- read_html(&quot;https://www.nytimes.com/interactive/2018/11/06/us/elections/results-minnesota-elections.html&quot;)
tbls &lt;- html_nodes(webpage, &quot;table&quot;)
# Shows all tables where Walz is matched
tbls %>% grep(&quot;Walz&quot;,.) %>% print() %>% tbls[.]</pre>

<p>This output two different tables that looked similar to the following…</p>

<pre>[1] 10 11
{xml_nodeset (2)}
[1] <table>...
[2] <table>...

<p>We want the second table with county results, which has an actual index of 11 as seen at the top of the output. This might not be the best way, but it is quick. Without the <code>%>% print() %>%</code> we just see [1] or [2], which would give us the wrong table.</p>

<pre>governer &lt;- html_table(tbls[11],fill = T)[[1]]
governer$County %&lt;>% tolower %>% gsub(&quot;\\.&quot;,&quot;&quot;,.)
governer$Walz %&lt;>% as.character %&lt;>% gsub(&quot;\\,&quot;,&quot;&quot;,.) %>% as.numeric
governer$Johnson %&lt;>% as.character %&lt;>% gsub(&quot;\\,&quot;,&quot;&quot;,.) %>% as.numeric
names(governer)[1] &lt;- &quot;subregion&quot;

# Merge with mn_county
mn_county &lt;- merge(mn_county,governer,by=&quot;subregion&quot;,all.x=T)
# Set the Margins
mn_county$margin &lt;- (mn_county$Walz / mn_county$Johnson) %>% -1</pre>

<p>The above was just to get the margins to diverge from 0 (which will be used later for <code>ggplot</code>.</p>

<pre>> min(mn_county$margin)
[1] -0.5754299
> max(mn_county$margin)
[1] 1.676873</pre>

<p>If you were wondering, that ~1.68 came from Ramsey County, voting 170,391 to 63,653 for Walz. The <code>min()</code> came from Morrison with 9,711 to 4,123 for Johnson.</p>
<p>Time to set the limits. I chose .3 in either direction, which seemed good.</p>

<pre>mn_county$margin[which(mn_county$margin>.3)] &lt;- .3
mn_county$margin[which(mn_county$margin&lt;(-0.3))] &lt;- (-0.3)</pre>

<h2>3D Rendering</h2>
<p>First up, the packages unique to this portion</p>

<pre>library(rayshader)
library(png)
library(raster)</pre>

<p>If you haven’t installed <code>rayshader</code>, use</p>

<pre>devtools::install_github(&quot;tylermorganwall/rayshader&quot;)</pre>

<h3>Generating Elevation</h3>
<p>To simplify the process of generating the elevation and color overlay I made two ggplots hiding typical elements. The first is for elevation.</p>

<pre>mn_3d &lt;- ggplot(data = mn_df, mapping = aes(x = long, y = lat, group = group)) + 
  geom_polygon(color = &quot;black&quot;, fill = &quot;black&quot;)</pre>

<p>This gives us just the rough outline for the state of Minnesota</p>
<p><img alt="MN ggplot" title="MN ggplot" src="https://i0.wp.com/rlang.io/wp-content/uploads/2018/11/MN-ggplot.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Plot the elevation using the county density.</p>

<pre>mn_3d &lt;- mn_3d + theme_nothing() + theme(plot.background = element_rect(fill = &quot;black&quot;)) +
  geom_polygon(data = mn_county, aes(fill = density), color = &quot;black&quot;) + 
  scale_fill_continuous(low = &quot;#010101&quot;,high = &quot;white&quot;) +
  theme(legend.position = &quot;none&quot;,
        axis.ticks = element_blank(),
        panel.grid = element_blank(),
        axis.text = element_blank()) +
  geom_polygon(color = &quot;black&quot;, fill = NA)  +
  labs(fill = &quot;&quot;) +  theme(plot.background = element_rect(fill = &quot;black&quot;))
mn_3d</pre>

<p><img alt="Elevation 2D" title="Elevation 2D" src="https://i1.wp.com/rlang.io/wp-content/uploads/2018/11/elevation-2d.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Save the plot, make sure to use the same dimensions for both elevation and color overlay.</p>

<pre>ggsave(filename = &quot;elevation-2d.png&quot;, plot = mn_3d, width = 6, height = 6)</pre>

<h3>Generating Color Overlay</h3>
<p>This is where the limits come into play.</p>

<pre>mn_gov &lt;- mn_3d + geom_polygon(data = mn_county, aes(fill = margin), color = &quot;black&quot;) + 
  scale_fill_continuous(limits = c(-.3,.3), low = &quot;#BD0000&quot;, high = &quot;#0040B8&quot;,
                        space = &quot;Lab&quot;, na.value = &quot;grey50&quot;, guide = &quot;colourbar&quot;) +
  theme(plot.background = element_rect(fill = &quot;#FFFFFF&quot;))
mn_gov
ggsave(filename = &quot;MN-Governor.png&quot;,plot = mn_gov,width = 6, height = 6)</pre>

<p><img alt="MN Governor" title="MN Governor" src="https://i2.wp.com/rlang.io/wp-content/uploads/2018/11/MN-Governor.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h3>Convert Plots to Matrix Values</h3>

<pre>raster::raster(&quot;elevation-2d.png&quot;) -> localtif
# And convert it to a matrix:
elmat &lt;- matrix(raster::extract(localtif,raster::extent(localtif),buffer=1000),
                nrow=ncol(localtif),ncol=nrow(localtif))</pre>

<p>For some reason ggplot seems to keep a bit of a border when using ggsave, so we will need to trim off a little from the matrix.</p>

<pre>> dim(elmat)
[1] 1800 1800</pre>

<p>I chose to trim by 10 on every side, for both elevation and color overlay.</p>

<pre>elmat &lt;- elmat[11:(nrow(elmat)-10),11:(ncol(elmat)-10)]</pre>

<p>Now we get this</p>

<pre>dim(elmat)
[1] 1780 1780</pre>

<h3>Do the Same for Color Overlay</h3>

<pre>ecolor &lt;- readPNG(&quot;MN-Governor.png&quot;)
ecolor &lt;- ecolor[11:(nrow(ecolor)-10),11:(ncol(ecolor)-10),1:4]
# Set the alpha value on 4th dimension (RGBA)
ecolor[,,4] &lt;- .9</pre>

<h2>The Payoff</h2>
<p>Here it is, the final piece from Rayshader. If you just want the elevation, remove the line for <code>add_overlay()</code>.</p>

<pre>elmat %>%
  sphere_shade(progbar = FALSE,texture = &quot;bw&quot;) %>%
  add_overlay(overlay = ecolor) %>%
  add_shadow(ray_shade(elmat,zscale=4000,maxsearch = 300,progbar = FALSE),0.7) %>%
  plot_3d(elmat, fov=30, theta=45, phi=25, windowsize=c(1024,1024), zoom=0.4,
          water=FALSE, waterdepth = 10, waterlinecolor = &quot;white&quot;, waterlinealpha = 0.5,
          wateralpha = 0.8,watercolor = &quot;lightblue&quot;)</pre>

<p>Rendering the water by changing that flag can be a good way to segment or filter out low density counties.</p>
<p><img alt="Election Results with Water" title="Election Results with Water" src="https://i0.wp.com/rlang.io/wp-content/uploads/2018/11/MN-Election-Results-Water.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h2>Some Useful Commands</h2>

<pre>library(rgl)
# Render Viewport to File
rgl.snapshot(&quot;MN-Election-Results-Water.png&quot;, fmt = &quot;png&quot;, top = F)
# Render Viewport with Simulated Depth of Field to Plot Viwer
render_depth(focus = .5, focallength = 85, fstop = 4)
# Render GIF (Can Take a While)
movie3d(spin3d(axis = c(0, 1, 0), rpm = 4), duration = 15, dir = getwd(), movie = &quot;render&quot;)</pre>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://rlang.io/create-3d-county-maps-using-density-as-z-axis/"> RLang.io | R Language Programming</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/TpVFwlk3h5A" height="1"></a></div></table></table></pre></table></pre>