For analyzing the effect of pollution, Mexico City was divided into five zones. The function get_zone_imeca can download data for each zone as measured in IMECAs as opposed to the original units that get_station_data uses.

## Auto-install required R packages
if (!require("pacman")) install.packages("pacman")
## Loading required package: pacman
pacman::p_load(aire.zmvm, dplyr, ggplot2)
## The mxmaps packages has to be installed from GitHub
## if (!require("devtools")) {
##     install.packages("devtools")
## }
## devtools::install_github("diegovalle/mxmaps")
library(mxmaps)

df <- zones
df$value <- df$zone

p <- mxmunicipio_choropleth(df, num_colors = 6,
                       zoom = df$region,
                       title = "Geographic Zones of Mexico City") +
  scale_fill_manual("Zone", values = c("Centro" = "#e41a1c", "Noreste" = "#377eb8",
                               "Noroeste" = "#4daf4a",
                               "Sureste" = "#984ea3",
                               "Suroeste" = "#ff7f00"))
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
plot(p)

Note that the standards for measuring PM10 and O3 in IMECAS changed in October 2014 and the function prints a message

# Download pm10 data since 2008 for all available zones ("TZ")
pm_10 <- get_zone_imeca(criterion = "MAXIMOS", # Can be MAXIMOS (daily maximum) or 
                                              # HORARIOS (hourly average)
                       pollutant = "PM10", # "SO2", "CO", "NO2", "O3", "PM10", 
                                           # "TC" (All pollutants)
                       zone = "TZ", # "NO", "NE", "CE", "SO", "SE", "TZ" (All zones)
                       start_date = "2010-01-01", # Can't be earlier than 2008-01-01
                       end_date = "2017-01-15") # Can be up to the current date
## Starting October 28, 2014 the IMECA values for O3 and PM10 are computed using NOM-020-SSA1-2014 and NOM-025-SSA1-2014
knitr::kable(head(pm_10))
date zone pollutant unit value
2010-01-01 NO PM10 IMECA 89
2010-01-02 NO PM10 IMECA 57
2010-01-03 NO PM10 IMECA 44
2010-01-04 NO PM10 IMECA 23
2010-01-05 NO PM10 IMECA 38
2010-01-06 NO PM10 IMECA 56

Plotting the data makes the change that took place in October 2014 really obvious:

# Plot the overall highest maximum pm10 level with trendline
ggplot(pm_10 %>% group_by(date) %>% summarise(max = max(value, na.rm = TRUE)),
                      aes(date, max, group = 1)) +
  geom_line(color = "darkgray", size = .2) +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 50), se = FALSE) +
  ggtitle("Daily maximum PM10 levels in IMECAS") +
  geom_vline(xintercept = as.numeric(as.Date("2014-10-28"))) +
  theme_bw()
## `summarise()` ungrouping output (override with `.groups` argument)