example_map_canada.Rmd

---
title: "Canada"
output: html_document
---

```{r setup, include=FALSE}
## Defaults for R chunks
knitr::opts_chunk$set(echo = TRUE, ## echo = TRUE means code will show
                      warning=FALSE, ## supress warnings and messages from R
                      message=FALSE,
                      # fig.path='Figs/', ## where to save figures
                      fig.width = 11)
## Add any R packages you require. 
## Here are some we will use in 811:
requires <- c("tidyverse", ## tidyverse includes dplyr and ggplot2
              "magrittr",
              "here",
              "sf",
              "geojsonio",
              "rgdal",
              "rmapshaper",
              "spdplyr")
## Install any you don't have
to_install <- c(requires %in% rownames(installed.packages()) == FALSE)
install.packages(c(requires[to_install], "NA"), repos = "https://cloud.r-project.org/" )
## Load all required R packages
library(tidyverse)
library(ggplot2); theme_set(theme_minimal())
library(magrittr)
library(here)
```

[From Kieran Healy's blog post](https://kieranhealy.org/blog/archives/2018/12/09/canada-map/).
My R Markdown file for this example is [here](https://github.com/judgelord/PS811/example_map_canada.Rmd).

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 tidyverse libraries provide. We start by loading the `sf` library, which will also bring a number of dependencies with it. We’ll also load a few packages that provide some functions we might use in the conversion and mapping process.

```{r}
library(geojsonio)
library(rmapshaper)
library(rgdal)
library(tidyverse)
library(spdplyr)
library(sf)
```

We make a function, theme_map(), 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.


```{r}
theme_map <- function(base_size=9, base_family="") {
    require(grid)
    theme_bw(base_size=base_size, base_family=base_family) %+replace%
        theme(axis.line=element_blank(),
              axis.text=element_blank(),
              axis.ticks=element_blank(),
              axis.title=element_blank(),
              panel.background=element_blank(),
              panel.border=element_blank(),
              panel.grid=element_blank(),
              panel.spacing=unit(0, "lines"),
              plot.background=element_blank(),
              legend.justification = c(0,0),
              legend.position = c(0,0)
              )
}

theme_set(theme_map())
```

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 (.shp) 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.

So, get the data from Statistics Canada. We’re going to use this link: http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2011-eng.cfm. From the linked page, choose as your options ArcGIS .shp file, then—for example—Census divisions and cartographic boundaries. 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 rgdal’s workhorse functions, readOGR.


```{r}
canada_raw <- readOGR(dsn = "data/gcd_000b11a_e", 
                      layer = "gcd_000b11a_e",
                      use_iconv=TRUE, 
                      encoding="CP1250")
```

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:


```{r}
canada_raw_json <- geojson_json(canada_raw)
canada_raw_sim <- ms_simplify(canada_raw_json)
```

Then we save the resulting GeoJSON file, which we can now work with directly from here on out:

```{r}
geojson_write(canada_raw_sim, file = here("data/canada_cd_sim.geojson"))
```

We only need to do the importing, converting, and thinning once. If you already have a GeoJSON file, you can just start here.

Read the GeoJSON file back in as an sf object. The .geojson file is included in the repository, so you can load the libraries listed above and then start from here if you want.

```{r}
canada_cd <- st_read(here("data/canada_cd_sim.geojson"), quiet = TRUE)
```

```{r canada}
canada_cd %>% 
ggplot() + 
  geom_sf(mapping = aes(fill = PRNAME)) +
  guides(fill = FALSE) 



head(canada_cd)
```

Notice the proj4string there in the metadata. Let us change that to Canada’s favorite projection, the Lambert Conformal Conic projection. 


```{r}
canada_cd %<>% st_transform(crs = "+proj=lcc +lat_1=49 +lat_2=77 +lon_0=-91.52 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")

head(canada_cd)
```

Finally, and just because I don’t have any census-district-level data 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.


```{r}
map_colors <-  RColorBrewer::brewer.pal(8, "Pastel1")
map_colors <- rep(map_colors, 37)
```

Instead, we’ll just map the fill to PRUID, i.e. the province level. But try mapping fill to CDUID instead (the census district ID), and see what happens.


```{r canada-projection}
canada_cd %>% 
ggplot(mapping = aes(fill = PRUID)) + # provence-level unit ID
  geom_sf( size = 0.1) + 
  scale_fill_manual(values = map_colors) + 
  guides(fill = FALSE) +
  theme_map() +
  theme(panel.grid.major = element_line(color = "white"),
        legend.key = element_rect(color = "gray40", size = 0.1))
```