In diegovalle/aire.zmvm: Download Mexico City Pollution, Wind, and Temperature Data

knitr::opts_chunk$set(fig.dim =  c(11,9),
                      fig.align = 'center',
                      collapse = TRUE,
                      comment = "#>",
                      out.width = '100%'
)

This document introduces you to tools for downloading data from each of the pollution, wind and temperature measuring stations in the Zona Metropolitana del Valle de México (greater Mexico City).

Core Functions

The basic set of functions provided by the package:

get_station_data and get_station_month_data download pollution, wind and temperature data for each of the measuring stations in the original units (ppb, µg/m³, etc).
get_station_imeca download pollution values for each station in IMECAs
get_zone_imeca download pollution data in IMECAs for each of the 5 geographic zones of Mexico City
get_latest_imeca download the latest pollution hourly maximums for each of the measuring stations.
idw360 inverse distance weighting modified to work with degrees, useful for wind data

Because of limitations with the SEDEMA website http://www.aire.cdmx.gob.mx/ (where all data downloaded by the aire.zmvm comes from), options for each of the function aren't entirely consistent but hopefully the following table will make them clearer.

| function | date range | Units | Wind / Tmp / RH | Earliest Date | Pollutants | Includes All Stations | Criterion | |---------------------------|----------------|----------|------------------|------------|------------|-------|-------| | get_station_data | years | Original | Yes | 1986 | SO2, CO, NO2, O3, PM10, PM25, WSP, WDR, TMP, RH | Yes | hourly, daily maximum, daily minimum | | get_station_month_data | 1 month | Original | Yes | 2005‑01 | SO2, CO, NO2, O3, PM10, PM25, WSP, WDR, TMP, RH | Yes | hourly, daily maximum, daily minimum | | get_station_imeca | 1 day | IMECA | No | 2009‑01‑01 | SO2, CO, NO2, O3, PM10 | No | hourly | | get_zone_imeca | 1 or more days | IMECA | No | 2008‑01‑01 | SO2, CO, NO2, O3, PM10 | Only zones | hourly, daily maximum | | get_latest_imeca | 1 hour | IMECA | No | Latest only| Maximum value of SO2, CO, NO2, O3, PM10 | No |latest hourly |

Map of Pollution Levels

Pollution is monitored by a series of stations located throughout the Mexico City Metro Area. Since there are different pollutants to monitor, each measured in different units, and each with different specified averaging periods, a common scale called IMECA was created to make them comparable. If pollutants reach 151 IMECA points or more the Mexico City government can declare a phase I smog alert and order cars off the streets and place limits on polluting industries.

On May 4, 2018 there was a smog alert, but there wasn no doble hoy no circula because metereological conditions were forecast to improve the next day. We'll show just how easy it is to make a map of pollution levels during the alert with this package.

First we load the packages we need to run this vignette and download the pollution data:

library("aire.zmvm")
library("lubridate")
library("dplyr")
library("ggplot2")
library("ggmap")
library("sp")
library("gstat")

# Download hourly data for all pollutants during May 04, 2018 
# The maximum Ozone level was reached at 17:00 PM local time,
# 16:00 with no daylight savings time
ldf <- lapply(c("PM10", "O3", "SO2", "NO2", "CO"), 
              function(x) get_station_imeca(x, "2018-05-04") %>% 
                filter(hour == 16) )

Here's what the data we just downloaded looks like:

knitr::kable(head(ldf[[2]][order(-ldf[[2]]$value), ], 10), row.names = FALSE)

You may have noticed that the code comments mention the local time as being 1 hour ahead of that returned by the function. That's because the hourly data returned by get_station_imeca comes in a format that is not affected by daylight savings time—i.e. 6 hours behind UTC (Etc/GMT+6). If, for example, you ever want to analyze pollution levels at peak commuting times, you'll want to convert the time to the local timezone of America/Mexico_City. You can use the following code to do the conversion:

format(as.POSIXct("2017-05-04 16:00", tz = "Etc/GMT+6"), 
       tz = "America/Mexico_City", 
       usetz = TRUE)

The package also includes a data.frame named stations that includes the latitude and longitude of all pollution measuring stations.

knitr::kable(head(stations, 5))

We will use this data to create a grid covering Mexico City and extrapolate the pollution values to each cell in the grid by inverse distance weighting from the locations of the measuring stations.

create_grid <- function(station_vec, pixels = 1000) {
  df <- stations[stations$station_code %in% station_vec,]
  geog <- df[,c("lat", "lon")]
  coordinates(geog) <- ~lon+lat
  proj4string(geog) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
  ## Add a margin surrounding the stations
  geog@bbox <- bbox(geog) + c(-0.05, -0.05, 0.05, 0.05)
  geog.grd <- makegrid(geog, n = pixels)

  grd.pts <- SpatialPixels(SpatialPoints(geog.grd))
  as(grd.pts, "SpatialGrid")
}


# Sometimes stations are taken offline for maintenance, exclude those
# that weren't reporting data
reporting_stations <- unique(na.omit(do.call(rbind, ldf))$station_code)
# 50x50 grid covering Mexico City
grid <- create_grid(reporting_stations, 15000)

plot(grid, main = "Grid covering Mexico City and Locations of\nPollution Measuring Stations (red)",
     col = "#666666", lwd = .2)
# Plot the locations of the stations
geog <- stations[stations$station_code %in% reporting_stations, c("lat", "lon")]
coordinates(geog) <- ~lon+lat
points(geog, pch = 22, col = "#333333", bg = "tomato1", cex = 1.2)

There are probably additive or interaction effects between pollutants, but to keep things simple we'll do like the Mexico City government when it reports pollution levels and simply ignore them. We'll take the maximum pollution value in each cell of the grid to create the map.

heatmap <- function(df, grid){
  if(nrow(df) == 0){
    return(data.frame(var1.pred = NA, var1.var = NA, lon = NA, lat = NA))
  }

  df <- left_join(df, stations, by = "station_code")
  df <- df[!is.na(df$value),]
  df <- df[,c("lat", "lon", "value")]
  coordinates(df) <- ~lon+lat
  # For radiation pollution the exponent should be 2
  # See http://www.sciencedirect.com/science/article/pii/S009830041200372X
  df.idw <- idw(value ~ 1, df, grid, idp = 2, debug.level = 0)

  idw = as.data.frame(df.idw)
  names(idw) <- c("var1.pred", "var1.var", "lon", "lat")

  idw
}


# Create the IDW heatmap for each pollutant
idw_tiles <- lapply(ldf, heatmap, grid)
# Calculate the maximum value across all pollutants
maxes <- apply(do.call(cbind, lapply(idw_tiles, 
                                     function(x) x[["var1.pred"]])), 
               1, max)
idw_tiles_max <- data.frame(var1.pred = maxes, 
                       lat = idw_tiles[[1]]$lat, 
                       lon = idw_tiles[[1]]$lon)

Plot

Now we will use the package ggmap to plot the locations of the stations and the pollution grid on top of a map. Note that you may have to install a beta version of ggmap by running: devtools::install_github("dkahle/ggmap", ref = "tidyup")

qmplot(x1, x2, data = data.frame(grid), geom = "blank", 
       maptype = "toner-lite", source = "stamen", zoom = 10)  +
  geom_tile(data = idw_tiles_max, aes(x = lon, y = lat, fill = var1.pred), 
            alpha = .8) + 
  scale_fill_gradient2("IMECA", 
                       low = "#abd9e9",
                       mid = "#ffffbf",
                       high = "#f46d43",
                       midpoint = 100) +
  geom_point(data = stations[stations$station_code %in% reporting_stations, ],
             aes(lon, lat), shape = 22, size = 1.6) +
  ggtitle(paste0("Pollution in Mexico City during the pollution",
                 " emergency of May 4, 2018"))