<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/canada-map/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://kieranhealy.org/blog/archives/2018/12/09/canada-map/"> R on kieranhealy.org</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>I taught my <a rel="nofollow" href="https://statisticalhorizons.com/seminars/public-seminars/data-visualization-fall18" target="_blank">Data Visualization seminar</a> in Philadelphia this past Friday and Saturday. It covers most of the content of <a rel="nofollow" href="http://socviz.co/" target="_blank">my book</a>, including a unit on making maps. The examples in the book are from the United States. But what about other places? Two of the participants were from Canada, and so here’s an example that walks through the process of grabbing a shapefile and converting it to a simple-features object for use in R. A self-contained R project with this code is <a rel="nofollow" href="https://github.com/kjhealy/canmap" target="_blank">available on GitHub</a>.</p>
<p>Mapping data (and, more importantly, converting map data into a format that we can tidily use) requires a bit more heavy lifting behind the scenes than the core <code>tidyverse</code> libraries provide. We start by (if necessary) installing the <code>sf</code> library, which will also bring a number of dependencies with it. We’ll also install a few packages that provide some functions we might use in the conversion and mapping process.</p>
<div>
<pre><code data-lang="r">

install.packages<span>(</span><span>&quot;sf&quot;</span><span>)</span>
install.packages<span>(</span><span>&quot;rgdal&quot;</span><span>)</span>
install.packages<span>(</span><span>&quot;geojsonio&quot;</span><span>)</span>
install.packages<span>(</span><span>&quot;spdplyr&quot;</span><span>)</span>
install.packages<span>(</span><span>&quot;rmapshaper&quot;</span><span>)</span>

</code></pre>
</div>
<p>Then we load the packages, along with others we’ve already installed (as part of the seminar).</p>
<div>
<pre><code data-lang="r">

<span>library</span><span>(</span>geojsonio<span>)</span>
<span>library</span><span>(</span>rmapshaper<span>)</span>
<span>library</span><span>(</span>rgdal<span>)</span>
<span>library</span><span>(</span>tidyverse<span>)</span>
<span>library</span><span>(</span>spdplyr<span>)</span>
<span>library</span><span>(</span>sf<span>)</span>
<span>library</span><span>(</span>socviz<span>)</span>


</code></pre>
</div>
<p>We make a function, <code>theme_map()</code>, that will be our ggplot theme. It consists mostly in turning off pieces of the plot (like axis text and so on) that we don’t need.</p>
<div>
<pre><code data-lang="r">

theme_map <span>&lt;-</span> <span>function</span><span>(</span>base_size<span>=</span><span>9</span><span>,</span> base_family<span>=</span><span>&quot;&quot;</span><span>)</span> <span>{</span>
    <span>require</span><span>(</span>grid<span>)</span>
    theme_bw<span>(</span>base_size<span>=</span>base_size<span>,</span> base_family<span>=</span>base_family<span>)</span> <span>%+replace%</span>
        theme<span>(</span>axis.line<span>=</span>element_blank<span>(),</span>
              axis.text<span>=</span>element_blank<span>(),</span>
              axis.ticks<span>=</span>element_blank<span>(),</span>
              axis.title<span>=</span>element_blank<span>(),</span>
              panel.background<span>=</span>element_blank<span>(),</span>
              panel.border<span>=</span>element_blank<span>(),</span>
              panel.grid<span>=</span>element_blank<span>(),</span>
              panel.spacing<span>=</span>unit<span>(</span><span>0</span><span>,</span> <span>&quot;lines&quot;</span><span>),</span>
              plot.background<span>=</span>element_blank<span>(),</span>
              legend.justification <span>=</span> <span>c</span><span>(</span><span>0</span><span>,</span><span>0</span><span>),</span>
              legend.position <span>=</span> <span>c</span><span>(</span><span>0</span><span>,</span><span>0</span><span>)</span>
              <span>)</span>
<span>}</span>

theme_set<span>(</span>theme_map<span>())</span>

</code></pre>
</div>
<p>Next we need to get the actual map data. This is often the hardest bit. In the case of Canada, their central statistics agency provides map shape files that we can use. The Shape File format (<code>.shp</code>) is the most widely-used standard for maps. We’ll need to grab the files we want and then convert them to a format the tidyverse can use. You won’t be able to run the code in the next few chunks unless the files are downloaded where the code expects them to be.</p>
<p>So, get the data from Statistics Canada. We’re going to use this link: <a rel="nofollow" href="http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2011-eng.cfm" target="_blank">http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2011-eng.cfm</a>. From the linked page, choose as your options <em>ArcGIS .shp file</em>, then—for example—<em>Census divisions</em> and <em>cartographic boundaries</em>. You’ll then download a zip file. Expand this zip file into the directory in your working folder named ‘data’. Then we can import the shapefile using one of <code>rgdal</code>’s workhorse functions, <code>readOGR</code>.</p>
<div>
<pre><code data-lang="r">

canada_raw <span>&lt;-</span> readOGR<span>(</span>dsn <span>=</span> <span>&quot;data/gcd_000b11a_e&quot;</span><span>,</span> layer <span>=</span> <span>&quot;gcd_000b11a_e&quot;</span><span>,</span>
                      use_iconv<span>=</span><span>TRUE</span><span>,</span> encoding<span>=</span><span>&quot;CP1250&quot;</span><span>)</span>
</code></pre>
</div>
<p>Note the options here. This information is contained in the documentation for the shape files. Now we’re going to convert this object to GeoJSON format and simplify the polygons. These steps will take a little while to run:</p>
<div>
<pre><code data-lang="r">

canada_raw_json <span>&lt;-</span> geojson_json<span>(</span>canada_raw<span>)</span>
canada_raw_sim <span>&lt;-</span> ms_simplify<span>(</span>canada_raw_json<span>)</span>

</code></pre>
</div>
<p>Then we save the resulting GeoJSON file, which we can now work with directly from here on out:</p>
<div>
<pre><code data-lang="r">
geojson_write<span>(</span>canada_raw_sim<span>,</span> file <span>=</span> <span>&quot;data/canada_cd_sim.geojson&quot;</span><span>)</span>
</code></pre>
</div>
<p>We only need to do the importing, converting, and thinning once. If you already have a GeoJSON file, you can just start here.</p>
<p>Read the GeoJSON file back in as an <code>sf</code> object. The <code>.geojson</code> file is included in the repository, so you can load the libraries listed above and then start from here if you want.</p>
<div>
<pre><code data-lang="r">
canada_cd <span>&lt;-</span> st_read<span>(</span><span>&quot;data/canada_cd_sim.geojson&quot;</span><span>,</span> quiet <span>=</span> <span>TRUE</span><span>)</span>

canada_cd

<span>## Simple feature collection with 293 features and 6 fields</span>
<span>## geometry type: MULTIPOLYGON</span>
<span>## dimension: XY</span>
<span>## bbox: xmin: -141.0181 ymin: 41.7297 xmax: -52.6194 ymax: 83.1355</span>
<span>## epsg (SRID): 4269</span>
<span>## proj4string: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs</span>
<span>## First 10 features:</span>
<span>## CDUID CDNAME CDTYPE PRUID PRNAME rmapshaperid geometry</span>
<span>## 1 4609 Division No. 9 CDR 46 Manitoba 0 MULTIPOLYGON (((-97.9474 50...</span>
<span>## 2 5901 East Kootenay RD 59 British Columbia / Colombie-Britannique 1 MULTIPOLYGON (((-114.573 49...</span>
<span>## 3 5933 Thompson-Nicola RD 59 British Columbia / Colombie-Britannique 2 MULTIPOLYGON (((-120.1425 5...</span>
<span>## 4 4816 Division No. 16 CDR 48 Alberta 3 MULTIPOLYGON (((-110 60, -1...</span>
<span>## 5 5919 Cowichan Valley RD 59 British Columbia / Colombie-Britannique 4 MULTIPOLYGON (((-123.658 48...</span>
<span>## 6 4621 Division No. 21 CDR 46 Manitoba 5 MULTIPOLYGON (((-99.0172 55...</span>
<span>## 7 4608 Division No. 8 CDR 46 Manitoba 6 MULTIPOLYGON (((-98.6436 50...</span>
<span>## 8 4811 Division No. 11 CDR 48 Alberta 7 MULTIPOLYGON (((-112.8438 5...</span>
<span>## 9 4802 Division No. 2 CDR 48 Alberta 8 MULTIPOLYGON (((-111.3881 5...</span>
<span>## 10 5951 Bulkley-Nechako RD 59 British Columbia / Colombie-Britannique 9 MULTIPOLYGON (((-124.4407 5...</span>

</code></pre>
</div>
<p>Notice the <code>proj4string</code> there in the metadata. We’re going to change that to Canada’s favorite projection, the Lambert Conformal Conic projection. See <a rel="nofollow" href="https://www.statcan.gc.ca/pub/92-195-x/2011001/other-autre/mapproj-projcarte/m-c-eng.htm" target="_blank">this discussion</a> for details.</p>
<div>
<pre><code data-lang="r">

canada_cd <span>&lt;-</span> st_transform<span>(</span>canada_cd<span>,</span> crs <span>=</span> <span>&quot;+proj=lcc +lat_1=49 +lat_2=77 +lon_0=-91.52 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs&quot;</span><span>)</span>

canada_cd

<span>## Simple feature collection with 293 features and 6 fields</span>
<span>## geometry type: MULTIPOLYGON</span>
<span>## dimension: XY</span>
<span>## bbox: xmin: -2529065 ymin: 5816874 xmax: 2793734 ymax: 10396190</span>
<span>## epsg (SRID): NA</span>
<span>## proj4string: +proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-91.52 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs</span>
<span>## First 10 features:</span>
<span>## CDUID CDNAME CDTYPE PRUID PRNAME rmapshaperid geometry</span>
<span>## 1 4609 Division No. 9 CDR 46 Manitoba 0 MULTIPOLYGON (((-457449.9 6...</span>
<span>## 2 5901 East Kootenay RD 59 British Columbia / Colombie-Britannique 1 MULTIPOLYGON (((-1628202 69...</span>
<span>## 3 5933 Thompson-Nicola RD 59 British Columbia / Colombie-Britannique 2 MULTIPOLYGON (((-1838098 74...</span>
<span>## 4 4816 Division No. 16 CDR 48 Alberta 3 MULTIPOLYGON (((-988280.5 7...</span>
<span>## 5 5919 Cowichan Valley RD 59 British Columbia / Colombie-Britannique 4 MULTIPOLYGON (((-2253700 71...</span>
<span>## 6 4621 Division No. 21 CDR 46 Manitoba 5 MULTIPOLYGON (((-461578.4 7...</span>
<span>## 7 4608 Division No. 8 CDR 46 Manitoba 6 MULTIPOLYGON (((-500377.6 6...</span>
<span>## 8 4811 Division No. 11 CDR 48 Alberta 7 MULTIPOLYGON (((-1350599 73...</span>
<span>## 9 4802 Division No. 2 CDR 48 Alberta 8 MULTIPOLYGON (((-1372235 69...</span>
<span>10</span>  <span>5951</span> Bulkley<span>-</span>Nechako     RD    <span>59</span> British Columbia <span>/</span> Colombie<span>-</span>Britannique            <span>9</span> <span>## MULTIPOLYGON (((-1921887 78...</span>

</code></pre>
</div>
<p>Finally, and just because I don’t have any census-district-level data to hand at the moment, we make a vector of repeated colors to fill in the map, for decoration only, if you want to color all the census divisions.</p>
<div>
<pre><code data-lang="r">

map_colors <span>&lt;-</span>  RColorBrewer<span>::</span>brewer.pal<span>(</span><span>8</span><span>,</span> <span>&quot;Pastel1&quot;</span><span>)</span>
map_colors <span>&lt;-</span> <span>rep</span><span>(</span>map_colors<span>,</span> <span>37</span><span>)</span>


</code></pre>
</div>
<p>Instead, we’ll just map the fill to <code>PRUID</code>, i.e. the province level. But try mapping <code>fill</code> to <code>CDUID</code> instead (the census district ID), and see what happens.</p>
<div>
<pre><code data-lang="r">

p <span>&lt;-</span> ggplot<span>(</span>data <span>=</span> canada_cd<span>,</span> 
            mapping <span>=</span> aes<span>(</span>fill <span>=</span> PRUID<span>))</span>
p_out <span>&lt;-</span> p <span>+</span> geom_sf<span>(</span>color <span>=</span> <span>&quot;gray60&quot;</span><span>,</span> 
                    size <span>=</span> <span>0.1</span><span>)</span> <span>+</span> 
  scale_fill_manual<span>(</span>values <span>=</span> map_colors<span>)</span> <span>+</span> 
  guides<span>(</span>fill <span>=</span> <span>FALSE</span><span>)</span> <span>+</span> 
  theme_map<span>()</span> <span>+</span> 
  theme<span>(</span>panel.grid.major <span>=</span> element_line<span>(</span>color <span>=</span> <span>&quot;white&quot;</span><span>),</span>
        legend.key <span>=</span> element_rect<span>(</span>color <span>=</span> <span>&quot;gray40&quot;</span><span>,</span> size <span>=</span> <span>0.1</span><span>))</span>

ggsave<span>(</span><span>&quot;figures/canada.pdf&quot;</span><span>,</span> p_out<span>,</span> height <span>=</span> <span>12</span><span>,</span> width <span>=</span> <span>15</span><span>)</span>

p_out

</code></pre>
</div>
<figure>
<img alt="Canada map" src="https://i0.wp.com/kieranhealy.org/files/misc/canada.png?w=456" data-recalc-dims="1"><figcaption>
<p>Map of Canada, colored by Province, showing 2011 Census District boundaries.</p>
</figcaption></figure>
<p>And there we have it. A map of Canada with census-division boundaries that you can join data to, in the same way as described in the “Draw Maps” chapter of <em><a rel="nofollow" href="http://socviz.co/" target="_blank">Data Visualization: A Practical Introduction</a></em>.</p>
	<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://kieranhealy.org/blog/archives/2018/12/09/canada-map/"> R on kieranhealy.org</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/lygFlSB6y38" height="1"></a></div>