Beijing PM 2.5 Concentration Line Chart

beijing_pm25.R

if (!require(pacman)){ install.packages("pacman") }
pacman::p_load(data.table, zoo, dygraphs)

# Data Source: http://www.stateair.net/web/historical/1/1.html
quality = rbind(
    fread("Beijing_2015_HourlyPM25_created20160201.csv", skip=3),
    fread("Beijing_2016_HourlyPM25_created20170201.csv", skip=3),
    fread("Beijing_2017_HourlyPM25_created20170705\ (3).csv", skip=3)
)

# There are 3 duplicate entries
sum(duplicated(quality[, get("Date\ (LST)")]))
# Use the latest entry
quality = quality[!duplicated(quality[, get("Date\ (LST)")], fromLast=T)]
    
# Impute missing values
quality[get("QC Name") != "Valid", Value := NA]
quality[, Value := na.locf(quality[, Value], na.rm=F)]

# Parse dates
quality[, date:=as.Date(strptime(quality[, get("Date (LST)")], "%m/%d/%Y %H:%M"), tz="Asia/Shanghai")]


# PM 2.5 24-hour Guideline: PDF Page 10
# http://apps.who.int/iris/bitstream/10665/69477/1/WHO_SDE_PHE_OEH_06.02_eng.pdf 

# US EPA 24-hour Guideline: 
# https://www.epa.gov/criteria-air-pollutants/naaqs-table

daily = quality[, .(median=as.numeric(median(Value)), avg=mean(Value),
                    upper=quantile(Value, probs=0.90), 
                    lower=quantile(Value, probs=0.1)), 
                by=date]


dygraph(daily[, .(date, lower, avg, upper)], 
        main="Beijing Daily PM2.5 Concentration") %>% 
    dyAxis("x", drawGrid=F) %>% dyAxis("y", drawGrid=F) %>%
    dySeries(c("lower", "avg", "upper"), label="PM2.5") %>%
    dyLimit(35, "US EPA 24-hr Guildline", labelLoc="right", color="green") %>%
    dyLimit(25, "WHO 24-hr Guildline", labelLoc="right", color="grey") %>%
    dyOptions(colors = RColorBrewer::brewer.pal(3, "Set1")) %>%
    dyRangeSelector(dateWindow=c(as.Date("2017-01-01"), max(daily$date))) %>% dyLegend(width = 200, show = "follow")