---
title: "Extract and Visualize the Population of an Urban City"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Extract and Visualize the Population of an Urban City}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
# The maps need rgeomex (locality boundaries) and several plotting packages.
# If any of them is missing, the code is shown but not run.
library(mapview)
mapviewOptions(basemaps = c("Esri.WorldGrayCanvas","Stadia.AlidadeSmoothDark","OpenStreetMap", "Esri.WorldImagery", "OpenTopoMap"))
needed <- c("rgeomex", "ggplot2", "cowplot", "ggspatial", "patchwork")
run <- all(vapply(needed, requireNamespace, logical(1), quietly = TRUE))
run_interactive <- run && all(vapply(c("mapview", "leafsync"), requireNamespace,
                                     logical(1), quietly = TRUE))
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = run,
  fig.width = 7,
  fig.height = 4
)
```

```{r setup}
library(popmex)
library(mapview)
```

The visualization of the population of an urban city is carried out in two
steps: first the population of the urban locality is extracted with
`extract_pop()`, and then it is mapped.

## Population of the metropolitan area of Guadalajara

### Extract the population

```{r guadalajara-extract}
gdl <- c("Guadalajara", "Tlaquepaque", "Zapopan", "Tonalá")

x <- popmex::extract_pop(year = 2023, cve_edo = "14", locality = gdl)
y <- popmex::extract_pop(year = 2022, cve_edo = "14", locality = gdl)

# total population in each year
c("2023" = sum(x$population), "2022" = sum(y$population))
```

### Static map

```{r guadalajara-static}
pop_map <- function(data, title, location = "bl") {
    ggplot2::ggplot() +
        ggplot2::geom_sf(data = data,
                         ggplot2::aes(fill = population),
                         col = "white",
                         lwd = 0.01) +
        ggplot2::scale_fill_viridis_c("Population") +
        ggplot2::labs(title = title) +
        cowplot::theme_map() +
        ggspatial::annotation_scale(style = "ticks", location = location)
}

library(patchwork)
pop_map(data = x, title = "2023") | pop_map(data = y, title = "2022")
```

### Interactive map

```{r guadalajara-interactive, eval = run_interactive}
a <- mapview::mapview(x, zcol = "population", layer.name = "Population 2023")
b <- mapview::mapview(y, zcol = "population", layer.name = "Population 2022")
leafsync::sync(a, b)
```

## Population of the metropolitan area of Mérida

### Extract the population

```{r merida-extract}
mid <- c("Merida", "Uman", "Kanasin", "Caucel", "Cholul")

x <- popmex::extract_pop(year = 2023, cve_edo = "31", locality = mid)
y <- popmex::extract_pop(year = 2022, cve_edo = "31", locality = mid)

c("2023" = sum(x$population), "2022" = sum(y$population))
```

### Static map

```{r merida-static}
pop_map(x, "2023", location = "br") | pop_map(y, "2022", location = "br")
```

### Interactive map

```{r merida-interactive, eval = run_interactive}
a <- mapview::mapview(x, zcol = "population", layer.name = "Population 2023")
b <- mapview::mapview(y, zcol = "population", layer.name = "Population 2022")
leafsync::sync(a, b)
```

## Population by municipality

The census datasets (`pop2000`, `pop2010`, `pop2020`) and the CONAPO
projections (`conapop`) give the population of each municipality, which is
the usual denominator of incidence rates.

```{r municipality, eval = TRUE}
# five most populated municipalities in 2020
head(pop2020[order(-pop2020$pobtot), c("CVEGEO", "pobtot", "pobfem", "pobmas")], 5)

# CONAPO projection for Mérida (31050) in 2026, both sexes
sum(conapop$pob_total[conapop$cvegeo == "31050" & conapop$ano == 2026])
```