---
title: "Maps and Geographical Data, Continued"
output:
  html_document:
    toc: yes
    code_folding: show
    code_download: true
---

<link rel="stylesheet" href="stat4580.css" type="text/css" />
<style type="text/css"> .remark-code { font-size: 85%; } </style>

```{r setup, include = FALSE}
source(here::here("setup.R"))
options(htmltools.dir.version = FALSE)
library(ggplot2)
knitr::opts_chunk$set(collapse = TRUE, fig.align = "center",
                      fig.height = 5, fig.width = 6)
library(lattice)
library(tidyverse)
theme_set(theme_minimal() +
          theme(text = element_text(size = 16)) +
          theme(panel.border = element_rect(color = "grey30", fill = NA)))
set.seed(12345)
```

```{r, include = FALSE}
source(here::here("maps2setup.R"))
```


## Simple Features

Using `geom_polygon` and `geom_point` works well for many simple
applications.

For more complex applications it is best to use special spatial data
structures.

An approach becoming widely adopted in many open-source frameworks for
working with spatial data is called _simple features_.

The package `sf` provides a range of tools for operating with simple
feature objects, such as:

* computing centroids;

* combining and intersecting objects;

* moving and expanding or shrinking objects.

Some useful vignettes in the `sf` package:

* [Simple Features for R](https://cran.r-project.org/package=sf/vignettes/sf1.html)

* [Reading, Writing and Converting Simple Features](https://cran.r-project.org/package=sf/vignettes/sf2.html)

* [Manipulating Simple Feature Geometries](https://cran.r-project.org/package=sf/vignettes/sf3.html)

Other references:

* Robin Lovelace, Jakub Nowosad, Jannes Muenchow (2019; 2nd Ed. 2025).
   _Geocomputation with R_, Chapman and Hall/CRC. Also available [on
   line](https://r.geocompx.org/).

* Mel Moreno and Mathieu Basille (2018). Drawing beautiful maps
  programmatically with R, sf and ggplot2, [Part 1:
  Basics](https://r-spatial.org/r/2018/10/25/ggplot2-sf.html),
  [Part 2:
  Layers](https://r-spatial.org/r/2018/10/25/ggplot2-sf-2.html),
  [Part 3:
  Layouts](https://r-spatial.org/r/2018/10/25/ggplot2-sf-3.html).

* Claudia A Engel (2019). [Using Spatial Data with
  R](https://cengel.github.io/R-spatial/).

* Andrew Heiss (2023). [Making Middle Earth maps with
  R](https://www.andrewheiss.com/blog/2023/04/26/middle-earth-mapping-sf-r-gis/).

* Taro Mieno (2023). [R as GIS for
  Economists](https://tmieno2.github.io/R-as-GIS-for-Economists/index.html).

A simple utility function for obtaining map data as sf objects:

<!-- ## nolint start: object_usage -->
```{r}
map_data_sf <- function(map, region = ".", exact = FALSE, crs = 4326, ...) {
    maps::map(map, region, exact = exact, plot = FALSE, fill = TRUE, ...) |>
        sf::st_as_sf(crs = crs) |>
        separate(ID, c("region", "subregion"), sep = ",") |>
        sf::st_transform(crs)
}
```
<!-- ## nolint end -->

Iowa county data:

```{r}
giowa_sf <- map_data_sf("county", "iowa")
```

This produces a special kind of data frame:

```{r}
class(giowa_sf)
```

A key part is the `geom` variable:

```{r}
head(giowa_sf, 4)
```

Johnson and Linn counties:

```{r}
filter(giowa_sf, subregion %in% c("johnson", "linn"))
```
```{r giowa-sf-jl, eval = FALSE}
filter(giowa_sf, subregion %in% c("johnson", "linn")) |>
    ggplot(aes(fill = subregion)) +
        geom_sf() +
        coord_sf()
```
```{r giowa-sf-jl, echo = FALSE, fig.width = 7}
```

A basic Iowa county map using `sf`:

```{r giowa-sf, eval = FALSE}
ggplot(giowa_sf) + geom_sf() + coord_sf()
```
```{r giowa-sf, echo = FALSE, fig.width = 7}
```

Adding wind turbine locations:

```{r wind-sf, eval = FALSE}
wt_IA_sf <-
    sf::st_as_sf(wt_IA,
                 coords = c("xlong", "ylat"),
                 crs = 4326)
ggplot(giowa_sf) +
    geom_sf(fill = "grey95") +
    geom_sf(aes(color = p_year),
            data = wt_IA_sf) +
    scale_color_viridis_c() +
    coord_sf()
```
```{r wind-sf, echo = FALSE, fig.width = 8}
```

It is possible to use `geom_point` but the `sf` coordinates may need
careful adjusting:

```{r wind-sf-point, eval = FALSE}
ggplot(giowa_sf) +
    geom_sf(fill = "grey95") +
    geom_point(aes(x = xlong, y = ylat,
                   color = p_year),
               data = wt_IA) +
    scale_color_viridis_c() +
    coord_sf() +
    labs(x = NULL, y = NULL)
```
```{r wind-sf-point, echo = FALSE, fig.width = 8}
```

With a Stamen/Stadia map background:

```{r wind-sf-stamen, eval = FALSE}
library(ggmap)
register_stadiamaps(StadiaKey, write = FALSE)
m <-
    get_stadiamap(c(-97.2, 40.4, -89.9, 43.6),
                  zoom = 8,
                  maptype = "stamen_toner")
ggmap(m) +
    geom_sf(data = giowa_sf,
            fill = NA,
            inherit.aes = FALSE) +
    geom_sf(aes(color = p_year),
            data = wt_IA_sf,
            inherit.aes = FALSE) +
    scale_color_viridis_c() +
    coord_sf() +
    labs(x = NULL, y = NULL)
```
```{r wind-sf-stamen, echo = FALSE, fig.width = 8, cache = TRUE, message = FALSE}
```

A choropleth map of population changes:

```{r pchanges-sf, eval = FALSE}
giowa_sf <- left_join(giowa_sf,
                      fipstab,
                      c("region", "subregion"))

giowa_pop_sf <-
    left_join(giowa_sf, cpop, "fips")

ggplot(giowa_pop_sf) +
    geom_sf(aes(fill = cpchange),
            show.legend = TRUE) +
    ggthemes::theme_map() +
    scale_fill_brewer(
        palette = "PRGn",
        drop = FALSE,
        guide = guide_legend(reverse = TRUE))
```
`show.legend = TRUE` is currently needed for `drop = FALSE` to work
properly.

```{r pchanges-sf, echo = FALSE, fig.width = 8}
```

Proportional symbol map for 2021 Iowa county population estimates:

```{r pop-syms-sf, eval = FALSE}
sf::sf_use_s2(FALSE)
iowa_pop_ctrds <-
    mutate(giowa_pop_sf,
           geom = sf::st_centroid(geom))
sf::sf_use_s2(TRUE)

ggplot(giowa_sf) +
    geom_sf() +
    geom_sf(aes(size = pop2021),
            data = iowa_pop_ctrds,
            show.legend = "point") +
    scale_size_area()
```

The use of `sf_use_s2()` can be avoided by using better map data.
```{r pop-syms-sf, echo = FALSE, warning = FALSE, message = FALSE, fig.width = 8}
```

A USA county map with a Bonne projection centered on Iowa City:

```{r pchanges-usa-sf, eval = FALSE}
gusa_sf <- map_data_sf("county")
gusa_sf <-
    left_join(gusa_sf,
              fipstab,
              c("region", "subregion"))
gusa_pop_sf <-
    left_join(gusa_sf, cpop, "fips")
ggplot(gusa_pop_sf) +
    geom_sf(aes(fill = cpchange), size = 0.1) +
    scale_fill_brewer(
        palette = "PRGn",
        na.value = "red",
        guide = guide_legend(reverse = TRUE)) +
    coord_sf(crs = "+proj=bonne +lat_1=41.6 +lon_0=-95")
```
```{r pchanges-usa-sf, echo = FALSE, fig.width = 8}
```

As an example of computation with `sf` objects we can revisit the
population contour plot.

First create an `sf` representation of the contours (there may be
better ways to do this):

<!--
based on
https://stackoverflow.com/questions/49526198/convert-a-fortified-data-frame-back-to-sf-object
-->

<!-- ## nolint start: object_usage -->
```{r}
ggpoly2sf <- function(poly,
                      coords = c("long", "lat"),
                      id = "group",
                      region = "region",
                      crs = 4326) {
    sf::st_as_sf(poly, coords = coords, crs = crs) |>
    group_by(!! as.name(id), !! as.name(region)) |>
    summarize(do_union = FALSE) |>
    sf::st_cast("POLYGON") |>
    ungroup() |>
    group_by(!! as.name(region)) |>
    summarize(do_union = FALSE) |>
    ungroup() |>
    sf::st_transform(crs)
}
```
<!-- ## nolint end -->

<!-- ## nolint start: object_usage -->
```{r}
sf_contourLines <- function(xyz, crs =  4326) {
    v <- contourLines(xyz)
    line_df <- function(i) {
        x <- v[[i]]
        x$group <- i
        as.data.frame(x)
    }
    lapply(seq_along(v), line_df) |>
        bind_rows() |>
        ggpoly2sf(coords = c("x", "y"),
                  region = "level",
                  crs = crs)
}
```
<!-- ## nolint end -->

```{r, message = FALSE}
dsi_sf <- sf_contourLines(dsi)
```

This reproduces our previous plot:

```{r contour-sf-1, eval = FALSE, fig.width = 7}
ggplot(giowa_sf) +
    geom_sf() +
    geom_sf(aes(fill = level),
            data = dsi_sf,
            alpha = 0.2)
```

```{r contour-sf-1, echo = FALSE, fig.width = 7}
```

With an outline of Iowa from the `maps` package we can use the
`sf::st_intersection()` function to show only the part of the contours
that is within Iowa.

```{r, message = FALSE, warning = FALSE}
iowa_sf <- map_data_sf("state", "iowa")
```

```{r, message = FALSE, warning = FALSE}
dsi_i_sf <-
    sf::st_intersection(iowa_sf, dsi_sf)
```

```{r conrour-sf-2, eval = FALSE}
ggplot(giowa_sf) +
    geom_sf() +
    geom_sf(aes(fill = level),
            data = dsi_i_sf,
            alpha = 0.2)
```

```{r conrour-sf-2, echo = FALSE, fig.width = 8, message = FALSE, warning = FALSE}
```

For the lower 48 states:

```{r, warning = FALSE, message = FALSE, fig.height = 4, class.source = "fold-hide"}
states_sf <- map_data_sf("state")
ds_sf <- sf_contourLines(ds)
ggplot(states_sf) +
    geom_sf() +
    geom_sf(aes(fill = level),
            data = ds_sf,
            alpha = 0.2)
```

```{r, warning = FALSE, message = FALSE, fig.height = 4, class.source = "fold-hide"}
sf::sf_use_s2(FALSE)
usa_sf <- sf::st_make_valid(map_data_sf("usa"))
ds_i_sf <- sf::st_intersection(usa_sf, ds_sf)
sf::sf_use_s2(TRUE)
ggplot(states_sf) +
    geom_sf() +
    geom_sf(aes(fill = level),
            data = ds_i_sf,
            alpha = 0.2)
```

The cleanup steps using `sf_use_s2()` and `st_make_valid()` can be
avoided by using better map data.


## Shape Files

The [ESRI](https://www.esri.com/en-us/home)
[_shapefile_ format](https://en.wikipedia.org/wiki/Shapefile) is a
widely used vector data format for geographic information systems.

Many GIS software systems support reading shapefiles.

Shapefiles are available from many sources; for example:

* The [TIGER](https://www.census.gov/programs-surveys/geography.html)
  site at the Census Bureau. The
  [`tigris`](https://cran.r-project.org/package=tigris) package
  provides an R interface.

* The [Natural Earth Project](https://www.naturalearthdata.com/).  The
  [rnaturalearth](https://cran.r-project.org/package=rnaturalearth)
  package provides an R interface.

Other formats and data bases are also available. Some references:

* Paul Murrell's _R Graphics, 2nd Ed._, Chapter 14.

* The [CRAN Task View](https://cran.r-project.org/web/views/) for
  [Analysis of Spatial Data](https://cran.r-project.org/view=Spatial)

* [Introduction to visualising spatial data in
  R](https://cran.r-project.org/doc/contrib/intro-spatial-rl.pdf).

* [Geocomputation with R](https://r.geocompx.org/).

<!--
Shape files for states and counties:

```{r, echo = FALSE, eval = FALSE}
options(tigris_use_cache = TRUE)
states <- tigris::states(resolution = "20m", cb = TRUE)
fp <- tibble::tibble(STATEFP = states[["STATEFP"]],
                     STNAME = states[["NAME"]])
cnty <- tigris::counties(resolution = "20m", cb = TRUE)
cnty <- left_join(cnty, fp, "STATEFP")
ggplot(filter(states, !(NAME %in% c("Hawaii", "Alaska", "Puerto Rico")))) +
    geom_sf()
ggplot(filter(cnty, !(STNAME %in% c("Hawaii", "Alaska", "Puerto Rico")))) +
    geom_sf()
```
-->

World map from Natural Earth:

```{r, fig.width = 10, class.source = "fold-hide"}
ne_zip <- "ne_10m_admin_0_countries.zip"
if (! file.exists(ne_zip)) {
    download.file(paste0("https://stat.uiowa.edu/~luke/data/", ne_zip),
                  ne_zip)
    unzip(ne_zip)
}
world <- sf::st_read("ne_10m_admin_0_countries.shp")
ggplot(world) +
    geom_sf(fill = "#E69F00") +
    theme(panel.background = element_rect(fill = "#56B4E950")) +
    coord_sf(crs = "+proj=moll")
```

State population isopleth map using shape files from the `tigris` package:

```{r, message = FALSE, results = FALSE}
options(tigris_use_cache = TRUE)
ssf <- tigris::states(resolution = "20m",
                      year = 2022,
                      cb = TRUE) |>
    filter(! NAME %in% c("Alaska", "Hawaii"))
```

```{r, message = FALSE}
ds_sf <- sf_contourLines(ds,
                         crs = sf::st_crs(ssf))
```

```{r state-isopleth-tigris, eval = FALSE}
ggplot(ssf) +
    geom_sf() +
    geom_sf(aes(fill = level),
            alpha = 0.2,
            data = sf::st_intersection(ssf,
                                       ds_sf))
```

```{r state-isopleth-tigris, echo = FALSE, warning = FALSE}
```


## Shape Files: Census Tracts

Census Tracts are small, relatively permanent subdivisions of a county.

Census tracts generally have  population sizes between 1,200 and 8,000 people.

Tracts for Iowa (state FIPS 19) for 2010:

<!-- include = FALSE to hide printing -->

```{r load-iowa-2010-tracts, include = FALSE}
tract_sf <- tigris::tracts(19, year = 2010)
```
```{r, eval = FALSE}
<<load-iowa-2010-tracts>>
```

```{r plot-iowa-cencus-tracts, eval = FALSE}
(p <- ggplot(tract_sf) + geom_sf())
```

```{r plot-iowa-cencus-tracts, echo = FALSE}
```

Johnson County:

```{r}
p %+% filter(tract_sf, COUNTYFP == "103")
```

Polk County:

```{r}
p %+% filter(tract_sf, COUNTYFP == "153")
```


## Proportion of Population Without Health Insurance

A lot of data is published at the census tract level.

One example is the proportion of the population that is uninsured
available from the [American Community
Survey](https://www.census.gov/programs-surveys/acs).

<!--
This example is adapted from two blog posts by Kevin Johnson
(unfortunately no longer available).

([part 1](http://www.kevjohnson.org/making-maps-in-r/),
[part 2](http://www.kevjohnson.org/making-maps-in-r-part-2/)).-->

Data from 2019 is available locally.

Download the data.

```{r, class.source = "fold-hide"}
dname <- "ACS_2019_5YR_S2701"
if (! file.exists(dname)) {
    data_url <- "https://stat.uiowa.edu/~luke/data"
    aname <- paste0(dname, ".tgz")
    download.file(paste(data_url, aname, sep = "/"), aname)
    untar(aname)
    file.remove(aname)
}
```

Read the data.
```{r, class.source = "fold-hide"}
acs_file <-
    file.path(dname,
              "ACSST5Y2019.S2701_data_with_overlays_2021-05-05T153143.csv")
acs_data_sf <- read.csv(acs_file, na.strings = "-")[-1, ] |>
    select(GEO_ID,
           name = NAME,
           population = S2701_C01_001E,
           uninsured = S2701_C05_001E) |>
    mutate(across(3:4, as.numeric)) |>
    mutate(GEO_ID = sub("1400000US", "", GEO_ID))
```

Merge and plot.

```{r plot-iowa-uninsured, eval = FALSE, class.source = "fold-hide"}
plot_data_sf <-
    left_join(tract_sf,
              acs_data_sf,
              c("GEOID10" = "GEO_ID"))

p <- ggplot(plot_data_sf) +
    geom_sf(aes(fill = uninsured)) +
    scale_fill_distiller(palette = "Reds",
                         direction = 1)
p
```
```{r plot-iowa-uninsured, echo = FALSE, fig.width = 9}
```

Johnson county:

```{r, class.source = "fold-hide"}
p %+%
    filter(plot_data_sf, COUNTYFP == "103")
```

Polk county:

```{r, class.source = "fold-hide"}
p %+%
    filter(plot_data_sf, COUNTYFP == "153")
```

For the state-wide map it may help to drop the census tract borders
and add county borders:

```{r plot-uninsured-2, eval = FALSE, class.source = "fold-hide"}
giowa_sf <-
    tigris::counties("iowa",
                     cb = TRUE,
                     resolution = "20m")
ggplot(plot_data_sf) +
    geom_sf(aes(fill = uninsured),
            color = NA) +
    geom_sf(fill = NA,
            color = "darkgrey",
            data = giowa_sf, size = 0.2) +
    scale_fill_distiller(palette = "Reds",
                         direction = 1) +
    coord_sf()
```

```{r plot-uninsured-2, echo = FALSE, fig.width = 9, message = FALSE}
```

The [`tidycensus`](https://walker-data.com/tidycensus/) package
provides a convenient way to work with census data, including the ACS.

```{r, include = FALSE}
keys <- as.data.frame(read.dcf(here::here("../APIKEYS")))
CensusKey <- keys$key[match("Census", keys$name)]
```
```{r}
library(tidycensus)
iowa_acs <- get_acs(key = CensusKey,
                    geography="tract",
                    variables = c("S2701_C01_001",
                                  "S2701_C05_001"),
                    state = "IA",
                    year = 2023,
                    geometry = TRUE,
                    output = "wide") |>
    select(GEOID,
           population = S2701_C01_001E,
           uninsured = S2701_C05_001E,
           geometry)

p <- ggplot(iowa_acs) +
    geom_sf(aes(fill = uninsured)) +
    scale_fill_distiller(palette = "Reds",
                         direction = 1)
```
```{r}
p %+% filter(iowa_acs, grepl("^19103", GEOID))
```

<!--
Might be nice to try mapgl:

https://walker-data.com/mapgl/

But mapbox seems to not play nicely with leaflet, plotly, and ggiraph.

```{r, include = FALSE}
keys <- as.data.frame(read.dcf(here::here("../APIKEYS")))
MapboxKey <- keys$key[match("CensusMapbox", keys$name)]
```
```{r, eval = FALSE}
m <- mapgl::mapboxgl(mapgl::mapbox_style("light"), iowa_acs, access_token = MapboxKey)
m
```
-->


## Interactive Maps: Leaflet

The [`leaflet` package](https://rstudio.github.io/leaflet/) can be
used to make interactive versions of these maps using the
[Leaflet](https://leafletjs.com/) JavaScript library for
interactive maps:

```{r uninsured-leaflet, fig.align = "center", class.source = "fold-hide"}
library(leaflet)
lf_sf <- leaflet(sf::st_transform(plot_data_sf, 4326), width = "70%")

addPolygons(lf_sf, fillColor = ~ colorNumeric("Reds", uninsured)(uninsured),
            weight = 0.5, color = "black", fillOpacity = 0.9, opacity = 1.0,
            highlightOptions = highlightOptions(color = "white",
                                                weight = 2,
                                                bringToFront = TRUE),
            label = ~ lapply(paste0(name, "<br/>",
                                    "pct: ", uninsured, "<br/>",
                                    "pop: ", population),
                             htmltools::HTML))
```


## Interactive Maps: Plotly

```{r uninsured-plotly, warning = FALSE, fig.height = 8, fig.width = 10, class.source = "fold-hide"}
p <- ggplot(plot_data_sf) +
    geom_sf(aes(fill = uninsured,
                text = paste0(name, "\n",
                              "pct: ", uninsured, "\n",
                              "pop: ", population)),
            color = NA) +
    geom_sf(fill = NA, color = "darkgrey", data = giowa_sf, size = 0.2) +
    scale_fill_distiller(palette = "Reds", direction = 1) +
    coord_sf() +
    ggthemes::theme_map()
plotly::ggplotly(p, tooltip = "text") |>
    plotly::style(hoverlabel = list(bgcolor = "white"))
```


## Interactive Maps: Ggiraph

```{r uninsured-ggiraph, fig.width = 10, class.source = "fold-hide"}
p <- ggplot(plot_data_sf) +
    ggiraph::geom_sf_interactive(
                 aes(fill = uninsured,
                     tooltip = paste0(name, "\n",
                                      "pct: ", uninsured, "\n",
                                      "pop: ", population)),
                 color = NA) +
    geom_sf(fill = NA, color = "darkgrey", data = giowa_sf, size = 0.2) +
    scale_fill_distiller(palette = "Reds", direction = 1) +
    coord_sf() +
    ggthemes::theme_map()
ggiraph::girafe(ggobj = p)
```


## Shape Files: Zip Code Boundaries

Data is often aggregated at a zip code level.

Shape files for zip code boundaries are available from the Census
Bureau TIGER project and can be accessed with the `zctas` function in
the `tigris` package.

<!-- year = 2020, the current default, doesn't seem to work at the moment -->
```{r zip_code_file_load, eval = FALSE, class.source = "fold-hide"}
library(sf)
library(tigris)
library(ggplot2)
z_sf <-
    zctas(cb = TRUE,
          starts_with = c("50", "51", "52"),
          class = "sf",
          year = 2020)
```

```{r, cache = TRUE, include = FALSE, ref.label = "zip_code_file_load"}
## separate chunk to avoid output
```

A map of the Iowa zip code regions:

```{r iowa-zip-code-map, eval = FALSE, class.source = "fold-hide"}
ggplot(z_sf) +
    geom_sf(fill = NA,
            color = "black",
            size = 0.1) +
    coord_sf()
```
```{r iowa-zip-code-map, echo = FALSE, fig.width = 9}
```

Zip codes within Johnson County:

```{r, echo = FALSE, eval = FALSE, class.source = "fold-hide"}
## the simpler form below seems to work
jc_sf <- map_data_sf("county", "iowa",
                     crs = st_crs(z_sf)) |>
    filter(subregion == "johnson") |>
    sf::st_transform(sf::st_crs(z_sf))
```
```{r, echo = FALSE, eval = FALSE}
##this works but maybe cleaner to use tigris::counties
jc_sf <- map_data_sf("county", "iowa") |>
    filter(subregion == "johnson") |>
    sf::st_transform(sf::st_crs(z_sf))
```
```{r, warning = FALSE}
jc_sf <- tigris::counties(19, resolution = "500k", cb = TRUE) |>
    filter(COUNTYFP == 103)
ggplot(sf::st_intersection(z_sf, jc_sf)) +
    geom_sf() +
    geom_sf_label(aes(label = GEOID20))
```

Zip codes that intersect Johnson County:

```{r, warning = FALSE, class.source = "fold-hide"}
filter(z_sf,
       sf::st_intersects(z_sf,
                         jc_sf,
                         sparse = FALSE)) |>
    ggplot() +
    geom_sf(color = "black", data = jc_sf) +
    geom_sf(fill = NA) +
    geom_sf_label(aes(label = GEOID20))
```


## A Tile Map of State Populations

One package available for constructing state-level tile maps is
`statebins`.

The package is based on `ggplot` and provides functions for handling
continuous and discrete data.

```{r, fig.width = 10, message = FALSE}
spops <- group_by(pep2021, state, year) |>
    summarize(pop = sum(pop)) |>
    ungroup() |>
    filter(year == 2021)

library(statebins)
ggplot(spops) +
    geom_statebins(aes(fill = pop, state = state)) +
    coord_equal() +
    theme_statebins() +
    scale_fill_distiller(palette = "Reds", direction = 1)
```


## A Hex Map of State Populations

Hexagon maps are an alternative.

```{r, fig.width = 10, class.source = "fold-hide"}
## download the shape file
hexlayer <- "us_states_hexgrid"
hexdir <- hexlayer
if (! dir.exists(hexdir)) {
    hexzip <- paste0(hexlayer, ".zip")
    data_url <- "https://stat.uiowa.edu/~luke/data"
    download.file(paste(data_url, hexzip, sep = "/"), hexzip)
    dir.create(hexdir)
    unzip(hexzip, exdir = hexdir)
    file.remove(hexzip)
}

## read shape file and merge with data
hex_pops <- sf::st_read(dsn = hexdir, layer = hexlayer) |>
    mutate(state = sub(" \\(.*", "", google_nam)) |>
    left_join(spops, "state")

ggplot(hex_pops, aes(fill = pop)) +
    geom_sf() +
    geom_sf_text(aes(label = iso3166_2,
                     color = ifelse(rank(pop) < 48, "dark", "light"))) +
    scale_color_manual(values = c(dark = "black", light = "white"),
                       guide = "none") +
    scale_fill_distiller(palette = "Reds", direction = 1) +
    ggthemes::theme_map() +
    theme(legend.position = "bottom")
```


## A Cartogram of State Populations

Cartograms are another approach.

Cartograms deform regions so their areas match values of a variable.

```{r, message = FALSE, fig.width = 8, class.source = "fold-hide"}
states_sf <- sf::st_as_sf(maps::map("state", plot = FALSE, fill = TRUE),
                          crs = 4326)
left_join(mutate(states_sf, state = as.character(ID)),
          mutate(spops,
                 state.abb = state.abb[match(state, state.name)],
                 state = tolower(state)),
          "state") |>
    sf::st_transform(3857) |>
    cartogram::cartogram_cont("pop") |>
    ggplot() +
    geom_sf() +
    geom_sf_text(aes(label = state.abb), size = 3) +
    ggthemes::theme_map()
```

Cartograms are sometimes used to show electoral college results.


## A Cartogram of 2016 Presidential Election Votes

Scaled by total number of votes in each state:

```{r, results = FALSE, class.source = "fold-hide"}
states_sf <- tigris::states(resolution = "20m", cb = TRUE)

spres2016 <- pres |>
    group_by(fp = fips_code %/% 1000) |>
    summarize(margin = sum(dem_2016 - gop_2016, na.rm = TRUE),
              total = sum(dem_2016 + gop_2016 + oth_2016, na.rm = TRUE)) |>
    ungroup() |>
    mutate(win = ifelse(margin > 0, "DEM", "GOP"))
```

```{r, message = FALSE, fig.width = 8, class.source = "fold-hide"}
mutate(states_sf, fp = as.numeric(STATEFP)) |>
    left_join(spres2016, "fp") |>
    left_join(data.frame(NAME = state.name, ABB = state.abb), "NAME") |>
    filter(! (NAME %in% c("Hawaii", "Alaska", "Puerto Rico"))) |>
    sf::st_transform(3857) |>
    cartogram::cartogram_cont("total") |>
    ggplot() +
    geom_sf(aes(fill = win)) +
    geom_sf_text(aes(label = ABB)) +
    scale_fill_manual(values = c(DEM = "blue", GOP = "red")) +
    ggthemes::theme_map()
```


## Linked Micromaps

The `micromap` package provides several functions for creating linked
micromaps.

```{r jc-micromap, eval = FALSE, class.source = "fold-hide"}
library(micromap)

maptab <- filter(tract_sf, COUNTYFP == 103) |>
    sf::as_Spatial() |>
    create_map_table(IDcolumn = "GEOID10")
acs_data_jc <- filter(acs_data_sf, grepl("Johnson", name))

lmplot(stat.data = acs_data_jc, map.data = maptab,
       grouping = 5, ord.by = "uninsured",
       panel.types = c("labels", "dot", "dot", "map"),
       panel.data = c("GEO_ID", "uninsured", "population", NA),
       map.link = c("GEO_ID", "ID"),
       panel.att = list(list(1, text.size = 0.6),
                        list(2, header = "Pct Uninsured"),
                        list(3, header = "Population")))
```

```{r jc-micromap, echo = FALSE, message = FALSE, fig.asp = 1.3, fig.width = 6.5}
```


## Odds and Ends: Including Alaska and Hawaii

Distance and size are challenging:

Mercator:
```{r}
ggplot(states_sf) + geom_sf()
```

Bonne at Iowa City latitude:
```{r}
ggplot(states_sf) +
    geom_sf() +
    coord_sf(crs = "+proj=bonne +lat_1=41.6 +lon_0=-95")
```

One of several options for obtaining maps that include Alaska and
Hawaii is the `urbnmapr` package available on GitHub.

```{r urbanmapr-states, message = FALSE, class.source = "fold-hide"}
## remotes::install_github("UrbanInstitute/urbnmapr")
library(urbnmapr)
states_sf <- get_urbn_map(map = "states",
                          sf = TRUE)
ggplot(states_sf) + geom_sf()
```

```{r urbanmapr-counties, message = FALSE, class.source = "fold-hide"}
counties_sf <- get_urbn_map(map = "counties",
                            sf = TRUE)
ggplot(counties_sf) + geom_sf()
```

2016 presidential election winners by county:

```{r urbanmapr-counties-pres, fig.width = 10, class.source = "fold-hide"}
## fix Oglala Lakota County SD FIPS
pres <- mutate(pres, fips = ifelse(fips_code == 46113, 46102, fips_code))

counties_sf <- mutate(counties_sf, fips = as.numeric(county_fips))
counties_sf <- left_join(counties_sf, pres, "fips")

NA_col <- "grey80"
ggplot(counties_sf) +
    geom_sf(aes(fill = win), size = 0.1) +
    scale_fill_manual(values = c(DEM = "blue", GOP = "red", "NA" = NA_col),
                      na.value = NA_col) +
    ggthemes::theme_map()
```

Another option is to use shape data and `shift_geometry()` from
the `tigris` package:

```{r tigris-shift, results = FALSE, class.source = "fold-hide"}
tcounties_sf <- tigris::counties(resolution = "20m", cb = TRUE) |>
    ## filter(STATE_NAME != "Puerto Rico") |>
    mutate(fips_code = as.numeric(paste0(STATEFP, COUNTYFP))) |>
    left_join(pres, "fips_code") |>
    tigris::shift_geometry()
```

```{r tigris-shift-plot, fig.width = 10, class.source = "fold-hide"}
ggplot(tcounties_sf) +
    geom_sf(aes(fill = win)) +
    scale_fill_manual(values = c(DEM = "blue", GOP = "red")) +
    ggthemes::theme_map()
```


## Odds and Ends: Dot Density Maps

Dot density map of 2012 and 2016 election results for Iowa counties (using
`tidycensus::as_dot_density()`):

```{r dot-density-election-iowa, cache = TRUE, class.source = "fold-hide"}
iowa_sf <- filter(counties_sf, state_abbv == "IA")
iowa_sf_long <- pivot_longer(iowa_sf,
                             contains("20"),
                             values_to = "count",
                             names_to = "vname") |>
    separate(vname, c("party", "year"), convert = TRUE) |>
    filter(party %in% c("dem", "gop")) |>
    pivot_wider(names_from = "party", values_from = "count") |>
    mutate(margin = abs(dem - gop), win = ifelse(dem > gop, "DEM", "GOP")) |>
    sf::st_as_sf()

## the GitHub version of tidycensus may be needed for as_dot_density
## remotes::install_github("walkerke/tidycensus")
library(tidycensus)
dd <- pivot_longer(iowa_sf_long, dem | gop,
                   values_to = "votes",
                   names_to = "party") |>
    as_dot_density("votes", values_per_dot = 1000, group = "party")

p <- ggplot(filter(iowa_sf_long, year >= 2012)) +
    geom_sf() +
    geom_sf(aes(color = party),
            data = filter(dd, year >= 2012),
            size = 0.2) +
    facet_wrap(~ year) +
    scale_color_manual(name = "Party", values = c(dem = "blue", gop = "red")) +
    coord_sf(crs = "+proj=bonne +lat_1=41.6 +lon_0=-95") +
    guides(color = guide_legend(override.aes = list(size = 5))) +
    ggthemes::theme_map() +
    theme(legend.position = "right")
```

```{r, echo = FALSE, fig.width = 14, fig.height = 6}
## needed to evaluate p if the previous chunk is cached
p
```

An alternative is to show the margin (DEM votes minus GOP Votes) as
proportional symbols:

```{r, message = FALSE, fig.width = 14, fig.height = 6, class.source = "fold-hide"}
ctr <- sf::st_centroid(iowa_sf_long)
ggplot(filter(iowa_sf_long, year >= 2012)) +
    geom_sf() +
    geom_sf(aes(size = margin, color = win),
            data = filter(ctr, year >= 2012)) +
    scale_color_manual(name = "Winner", values = c(DEM = "blue", GOP = "red")) +
    coord_sf(crs = "+proj=bonne +lat_1=41.6 +lon_0=-95") +
    facet_wrap(~year) +
    guides(color = guide_legend(override.aes = list(size = 5))) +
    ggthemes::theme_map() +
    theme(legend.position = "right")
```

Another option is choropleth maps:

```{r, fig.width = 14, fig.height = 6, class.source = "fold-hide"}

## margin/win
filter(iowa_sf_long, year >= 2012) |>
ggplot(aes(fill = win)) +
    geom_sf() +
    scale_fill_manual(name = "Winner", values = c(DEM = "blue", GOP = "red")) +
    coord_sf(crs = "+proj=bonne +lat_1=41.6 +lon_0=-95") +
    ggthemes::theme_map() +
    theme(legend.position = "right") +
    facet_wrap(~year)
```


## Odds and Ends: Package [`mapview`](https://r-spatial.github.io/mapview/)

```{r}
mapview::mapviewOptions(basemaps.color.shuffle = FALSE)
filter(ctr, year == 2012) |>
    mapview::mapview(zcol = "win", cex = "margin")
```

```{r}
filter(ctr, year == 2016) |>
    mapview::mapview(zcol = "win", cex = "margin")
```


## Odds and Ends: Accessing  US Governmet Data

A lot of data is made available via [DATA.GOV](https://data.gov/).

Data files are available in various formats.

Web APIs are also available for accessing many data sets.

R packages are available for using many of these APIs.

Some examples:

* [`blsAPI`](https://www.bls.gov/developers/home.htm) for the Bureau
  of Labor Statistics (BLS) API.
  

* [`censusapi`](https://github.com/hrecht/censusapi) and
  [`tidycensus`](https://walker-data.com/tidycensus/) for the Census
  Bureau API.
  

* [`rnoaa`](https://github.com/ropensci/rnoaa) for NOAA data.

A book on working with cencus data is
[available](https://walker-data.com/census-r/index.html).


## Reading

Chapter [_Draw maps_](https://socviz.co/maps.html#maps)
in [_Data Visualization_](https://socviz.co/).

Chapter [_Visualizing geospatial
  data_](https://clauswilke.com/dataviz/geospatial-data.html) in
  [_Fundamentals of Data
  Visualization_](https://clauswilke.com/dataviz/).

```{r, eval = FALSE, echo = FALSE}
## mapview of Iowa, Johnson County

## Dot density map for uninsured in Johnson County:
plot_data_jc <- filter(plot_data_sf, COUNTYFP10 == 103)

dd <- mutate(plot_data_jc,
             insured = population * (1 - uninsured / 100),
             uninsured = population * (uninsured / 100)) |>
    pivot_longer(uninsured | insured,
                 values_to = "count",
                 names_to = "status") |>
    as_dot_density("count", values_per_dot = 100, group = "status")

ggplot(plot_data_jc) +
    geom_sf() +
    geom_sf(aes(color = status, size = status), data = dd) +
    scale_size_manual(values = c(insured = 1, uninsured = 1.5))
```
