---
author: "Sean Laverty"
title: "Class 10 notes"
date: "`r Sys.Date()`"
format: html
editor: source
---
## Packages

For personal projects, I like to use this to ensure that moving between computers and such, I have the packages I need.  Unpacking it from the inside out is helpful for understanding how it works.  Essentially, if the package name I think I will need is not in the list of installed packages, I install it.

```{r}
#| label: install-libraries
#| echo: true
if(!("maps" %in% row.names(installed.packages()))){install.packages("maps")}
if(!("mapproj" %in% row.names(installed.packages()))){install.packages("mapproj")}
if(!("mapdata" %in% row.names(installed.packages()))){install.packages("mapdata")}
```

Now load the packages, for some reason, `mapdata` didn't like being listed with the others.
```{r}
#| label: load-libraries
#| echo: true
library(maps, mapproj)
library(mapdata)
library(ggplot2)
```

## Selecting available data

We can make a world map in @fig-world-map using the command `map()`.  Taking in no additional arguments, by defaults it will produce a dull world map.

```{r}
#| label: fig-world-map
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "This is a boring world map."
#| echo: false
map()
```

Setting `database = 'state'` will give a United States map of the lower 48 states, as in @fig-us-state-map.

```{r}
#| label: fig-us-state-map
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "This is a boring US state map."
#| echo: true
map(database = 'state')
```

Making the same type of change, this time to a value of `'county'`, we get a county map as in @fig-us-county-map.  For a little flourish, we first plot the county boundaries in gray, and overplot black state boundaries.  This shows the tremendous variability, on average, of county area as you move East to West.  There are 30 counties, mostly out West, with areas larger than the entire state of New Hampshire.  More on that with a few special commands later.
```{r}
#| label: fig-us-county-map
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "This is a thrilling US county map."
#| echo: true
map(database = 'county', col = "gray")
map('state', add = TRUE)
```
We can continue to "zoom in" by setting `region = 'oklahoma'`.  There are a few things to note about the syntax for the map shown in @fig-oklahoma-county-map, first, or `region` is actually the _state_ name, additionally it is spelled with all lowercase.

```{r}
#| label: fig-oklahoma-county-map
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "This is a familiar (?) Oklahoma county map."
#| echo: true
map(database = 'county', region = 'oklahoma')
```
## Color and creativity

After checking the help file for `map()`, we saw some background information about the `region` argument, in somewhat of an attempt to understand it, we tested a sample argumend involving the exclusion of Svalbard island from the map of Norway.  After first making the map, we added a `red` fill on a `skyblue` background.  To highlight the distinction between the two localities (ugh, "regions") we overplot mainland Norway ("excluding" Svalbard Island) by highlighting with `yellow`.  With a little work we should be able to add a suitable legend.

One sort of important idea here is that syntax to exclude Svalbard.  That is based on Perl regular expressions.  These can be used for **incredibly** powerful text editing/mining, but that would be the subject of nearly an entire day if not course.

```{r}
#| label: fig-norway-mainland
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A Map of Norway (in red) with Svalbard island highlighted (in yellow)."
#| echo: true
par(bg = "skyblue")
map(region ="Norway", col = "red", fill = TRUE)
map(region ="Norway(?!:Svalbard)", col = "yellow", fill = TRUE, add = TRUE)
```

Interestingly enough, the above approach doesn't seem to work for highlighting Chinese provinces.  A failure seen in @fig-china.
```{r}
#| label: fig-china
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A map of China."
#| echo: true
par(bg = "skyblue")
map(region ="China", col = "red", fill = TRUE)
map(region ="China(?!:Qinghai)", col = "yellow", fill = TRUE, add = TRUE)
```
## Maps as outlines

To perhaps better understand maps, we can look at the state outline data. To do so we will use the function `map_data()` from the `ggplot2` library.  It makes conveniently available the data for state outlines. A bit of investigation (creative use of tables, or simply reading the documentation) shows that the data shows latitude and longitude, `group` represents the spatial unit of a given state, or more generally its `region`. 

```{r}
#| label: fig-us-map-scatter
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A plot of the geographic coordinates of the states of the United States, showing that the input data is points that trace boundaries."
#| echo: true
usStates <- map_data("state")
plot(usStates[, c("long", "lat")], pch = 19, cex = 0.1, type = 'p')
```

Had we set the argument of `map_data()` to `"county"`, those names would fill `subregion` column.  Now each *county* is a `group`.  Not just disjoint lands like islands or detached peninsulas.  The names `region` and `subregion` actually seem like great argument names after some thought.

To understand some of the additional variables a bit better, on-the-fly, we made a plot of the outline of Colorado, a simple (it sure seems) rectangle. By rotating the entries of `order` by 45 degrees, we can read the numbers a bit better and confirm that the maps is essentially "connect-the-dots" on a state-by-state basis.
```{r}
#| label: fig-colorado-map-scatter
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A plot of the geographic coordinates of the states of Colorado, showing that the input data is points that trace boundaries."
#| echo: true
plot(usStates[usStates$region == "colorado", c("long", "lat")], pch = 19, cex = 0.5, type = 'p')
text(usStates[usStates$region == "colorado", "long"], usStates[usStates$region == "colorado", "lat"], usStates[usStates$region == "colorado", "order"], srt = 45)
```

A more interesting example might be a state like Washington, or North Carolina.  Again, make a state-by-state table of `group` numbers to see why. As @fig-washington-map-scatter confirms, the outline of Washington is created by five froups of points.  In this more interesting case, it seems more helpful to color points by `group` rather than label points by `order`, as we had done in class.

```{r}
#| label: fig-washington-map-scatter
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A plot of the geographic coordinates of the states of Colorado, showing that the input data is points that trace boundaries."
#| echo: true
cols <- usStates[usStates$region == "washington", "group"]; cols <- cols - min(cols) + 1
plot(usStates[usStates$region == "washington", c("long", "lat")], pch = 19, cex = 0.5, type = 'p', col = cols)
```

Given that data, we can use the code printed with @fig-us-gg-map-scatter to plot a map. Our approach is to plot the map as a collection of filled polygons.  Interestingly, the seemingly redundant `group = group` argument, plays a large role. If we omit it, states with more than one land mass will be drawn having those land masses connected in odd ways.
```{r}
#| label: fig-us-gg-map
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A superfluously colorful plot of the United States."
#| echo: true
ggplot(data = usStates, 
       mapping = aes(x = long, y = lat,  
       group = group, fill = region)) + 
       geom_polygon(color = "black", linewidth = 0.1) + 
       guides(fill = FALSE)
```
To be clear about the code, there is an argoment to `ggplot()` called `mapping`, this has nothing to do with the fact that we are map-making.  This refers to the mapping of data to various properties of the eventual graph (e.g., soze, color, style, shape, depending, in general, on the geometry). The effect of this small change, not plotting by groups within states, is interesting (procedurally). But, as @fig-wash-gg-map shows, it is also completely wrong. As the polygons intersect the filled regions reverse.  This is the more colorful effect of what we saw with @fig-us-map-scatter when we set `type = 'l'`.

```{r}
#| label: fig-wash-gg-map
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A superfluously colorful plot of the United States."
#| echo: true
ggplot(data = subset(usStates, region == "washington"), 
       mapping = aes(x = long, y = lat,  
       fill = "pink")) + 
       geom_polygon(color = "black", linewidth = 0.1) + 
       guides(fill = FALSE)
```

Having learned a bit more about syntax, we can do a bit to improve the appearance by turning to a more accepted map projection.  The function `coord_map()` with arguments as shown allows us to change to one of many projections, some with their own necessary arguments, like `lat0` and `lat1`.  These change the scaling a bit, but it would require some systematic experimentation or careful reading to get a feel for their roles.  They seem like acceptable choices based on existing books and work by experts in the field.

```{r}
#| label: fig-us-gg-map-albers
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A superfluously colorful plot of the United States, but with better projection."
#| echo: true
ggplot(data = usStates, 
       mapping = aes(x = long, y = lat,  
       group = group, fill = region)) + 
       geom_polygon(color = "black", linewidth = 0.1) +
       coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
       guides(fill = FALSE)
```
## Sensible map coloring

Now we can start the real work.  By importing a second set of data here describing the results of a survey on regional variation in drink terminology, we can color states using elements of the data. Read more about the survey^[here: <https://www.businessinsider.com/soda-pop-coke-map-2018-10>, here: <https://popvssoda.com/>, or here: <https://popvssoda.com/statistics/ALL.html>.]
```{r}
#| echo: true
dat <- read.delim("./data/drink-names.csv", header = TRUE, sep = '\t')
head(dat)
```

The datasets share the identifier of state name, called `"region"` in the state map data and `"Region"` in the survey data.  These datasets (or arbitrary columns of them) can be merged together by looking where the rows match in the name of the state.  While it would have been just as easy to open the file and change case on either one of those, this is a nice opportunity to illustrate the flexibility of `merge()`.  Since we are overwriting the original map data variable, just be careful not to execute this column more than onece in sequence or you will re-merge additional copies of the survey data to teh original data and accumulate appended `.x`, `.y`, or worse `.x.y` and so on.  The easy fix is then just to start over by redefining the original merged object a few steps back.  By merging in this direction we assign all survey values to each component of the state outline data. To match, we begin by converting the case of the actual entries, so that `Alabama` will match with `alabama`, and so on.  Strangely the variable _names_ as well as their _*entries*_ differed in case.

```{r}
#| echo: true
dat$Region <- tolower(dat$Region)
usStates <- merge(usStates, dat, by.x = "region", by.y = "Region")
head(usStates)
```
The fact that there is one line of survey values per state means that those values repeat for _each_ data point corresponding to a point on the boundary of each state -- that is over 1000 points for Texas, with a large size and complicated perimeter.  The duplication does not cause harm here, and it is an unavoidable step in specifying `aes()`.

Though this does not make the best or most interesting use of the data, we can color by the statewide percent of the total population who responded as survey participants -- this give _no_ information about their actual answer in terms of preference.  Ultimately we might care more about their answers, but this builds a good example. The percent participating colors each state, with an automatically generated color scheme.  By comparing the map in @fig-us-gg-map-albers-pop to data -- always a good idea, it's too easy to make mistakes, misinterpret a value in data, or miss out on an opportunity to learn something -- we can verify, in the absence of a legend, that darker colors are lower percents and brighter colors are higher.

```{r}
#| label: fig-us-gg-map-albers-pop
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A thoughtfully colorful plot of the United States, but with better projection."
#| echo: true
ggplot(data = usStates, 
  mapping = aes(x = long, y = lat,  
  group = group, fill = Percent)) + 
  geom_polygon(color = "gray90", linewidth = 0.1, show.legend = TRUE) + 
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) + 
  guides(fill = FALSE)
```

There were a lot of only slightly different colors in @fig-us-gg-map-albers-pop, which did not necessarily help.  We our readers could get distracted by too subtle of details, perhaps. For simplicity we will omit the first line of data, which actually correspnds to the marginal summary.  It is easy to recreated this with `apply(dat[, -1], 2, sum)`.  We would not be able to calculate the mean of the first column, words, so we omit it with the syntax `[, -1]`.  Visual analysis of the basic histogram, suggests, if only for the purposes of interpretation, we could look at binning the percent values at each even percent up to $5\%$ and then possibly having a "catch-all" category like "$>5\%$".

```{r}
#| label: fig-us-pop-hist
#| fig-width: 5
#| fig-height: 5
#| fig-cap: "A histogram of participation percent values for the 50 participating states."
#| echo: true
dat <- dat[-1, ]
hist(dat$Percent)
##head(dat)
```

It's worth experimenting with the layers of commands in the first line, but essentially `cut()` is a very slick way to decide which bin each value belongs to.  Essentially this converts a number into a category. The resultung variable, `bin`, can be used to color.  Explore with the bin spacings, perhaps every $2\%$ from $0\%$ to $10\%$, since in this data nothing exceeds $10%$ and that might seem like an intuituve cutoff. There is not a _rule_ for making this decision.  Look at data on the side, in a number of ways, avoid choices you would have a hard time explaining, and be honest about what you see.  

```{r}
#| label: fig-us-gg-map-albers-bin
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "A thoughtfully colorful plot of the United States, but with better projection."
#| echo: true
usStates$bin2 <- cut(usStates$Percent, c(0, 1, 2, 3, 4, 5, 100))
usStates$bin <- as.numeric(usStates$bin2)

ggplot(data = usStates, 
  mapping = aes(x = long, y = lat,  
  group = group, fill = bin)) + 
  geom_polygon(color = "gray90", linewidth = 0.1, show.legend = TRUE) + 
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  scale_fill_continuous(
#    labels = levels(usStates$bin2),
    labels = c("0-1%","1-2%", "2-3%", "3-4%", "4-5%", "5+%"),
    guide = guide_legend(title = "Percent participating")
  )
```

## The controversy

Ideally we would soon start looking at coloring using actual interesting values in the data.  But, to avoid the xkcd "all maps are just population maps" pitfall (change `fill = bin` to `fill = Total` above to get a sense of what the population map would look like.), we wouldn't want to use any of the raw vote _numbers_. We would want to convert them to something like "per 100,000" to standardize.  Even then there are issues with other kinds of data where, for good reasons, like privacy, certain data, like suicide incidences, are not even included in certain published reports because in small populations that information could become personally-identifiable.  One reasonable strategy would be to look at shares of the vote.  With transparency you could possibly layer maps with contrasting, but pleasantly and perceptibly overlapping colors or filled contours.  Or, variations on faceting or small-multiples with beverage-as-panl. 

Go experiment.  As you can see, making the maps alone was a task, now carefully superimposing or incorporating spatial data becomes the next step.
