~thrill_seeker/dataviz-challenge-2019

8939ff098768daec1487cb03f34c1ea76e9f41bc — ThrillSeeker 4 months ago a195c99
First complete version of .Rmd documents; pending changing verbs to present
3 files changed, 549 insertions(+), 28 deletions(-)

M analysis.R
A dataviz_cat.Rmd
R dataviz2019_eng.Rmd => dataviz_eng.Rmd
M analysis.R => analysis.R +0 -1
@@ 188,7 188,6 @@ predictions_dt <- daily_dt[, apply_models(dt=.SD, model_dates=apply_model_dates)
daily_dt <- rbind(daily_dt, predictions_dt)

# Data for app - select and save dataset

stations_dt <- unique(quality_dt[, .(estacio, nom_districte, latitud, longitud)])
daily_dt <- merge(daily_dt, stations_dt, by=c('estacio'))
daily_dt <- daily_dt[, estacio:=NULL]

A dataviz_cat.Rmd => dataviz_cat.Rmd +412 -0
@@ 0,0 1,412 @@
---
title: 'Anàlisi dels nivells de contaminació a Barcelona'
output:
  html_document
---

# Introducció

Aquest anàlisi va ser realitzat per a poder participar al [World Dataviz Challenge 2019](https://opendata-ajuntament.barcelona.cat/ca/data-viz-kobe-home).
Amb la recent preocupació envers temes medioambientals, i tenint en compte la recomanació explícita de la competició per a considerar aquest tema, l'elecció era prou clara.
[R](https://www.r-project.org/), un dels llenguatges de programació més estesos en l'anàlisi de dades (i el meu preferit), va ser l'escollit per a desenvolupar aquest anàlisi.

L'anàlisi seguí els següents passos:
1. Ingesta: Ingestió de dades sobre contaminació.
2. Preprocessament: Processar les dades fins que estiguin estructurades per a l'anàlisi.
3. Exploració: Extreure patrons i evolucions de la contaminació utilitzant gràfiques.
4. Modelització: Entrenar i aplicar un model per tal de predir futurs nivells de contaminació.
5. Aplicació: Desenvolupar una aplicació que mostra l'evolució del nivell de contaminació per estació.

El codi font per a l'anàlisi, documentació i aplicació es pot trobar a [sourcehut](https://git.sr.ht/~thrill_seeker/dataviz-challenge-2019).
L'aplicació es pot trobar al [portal de shiny](x).

# Ingesta

En primer lloc, era necessari carregar les dades. Aquestes van ser extretes dels següents enllaços:
- [Qualitat de l'aire](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/qualitat-aire-detall-bcn)
- [Estacions de mesura de contaminació](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/qualitat-aire-estacions-bcn)
- [Contaminants](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/contaminants-estacions-mesura-qualitat-aire )

Tot i que existeix la possibilitat d'utilitzar una API directament per a aconseguir les dades, en aquest cas vaig preferir descarregar directament els fitxers en format .csv.

Els següents paquets van ser carregats pel desenvolupament.

```{r package-loading}
library(data.table)
library(dplyr)
library(forecast)
library(ggplot2)
library(leaflet)
library(Metrics)
library(scales)
library(shiny)
library(shinydashboard)
library(shinyWidgets)
```

De la mateixa manera que s'especifica a la URL que conté les dades sobre qualitat de l'aire, la estructura dels datasets havia canviat a partir de l'Abril del 2019.
Vaig decidir definir dues parts diferents d'ingesta: 
1. La *vella* que llegeix les dades a partir de juny de 2018 (primer arxiu disponible) fins a març de 2019 (l'últim arxiu amb la mateixa estructura).
2. La *nova* que llegeix les dades a partir d'abril de 2019 fins a setembre del 2019 (l'últim arxiu pel moment). Variables auxiliars van ser definides també.

```{r variable-loading}
# Define variables
input_path <- 'data/input'
output_path <- 'data/output'
filename <- '*_qualitat_aire_bcn.csv'
limit_dates <- list(min=as.Date('2018-06-01'), max=as.Date('2019-09-30'))

# Get files, but take into account that from April '19 the structure is different
files_to_read <- list.files(path=input_path, pattern=filename, ignore.case=TRUE)
filter <- as.numeric(gsub(pattern='[^0-9]', replacement='', x=files_to_read))
old_files <- files_to_read[!(filter >= 201904)]
new_files <- files_to_read[(filter >= 201904)]
```

Then it was possible to load the data and inspect both structures.
For most of data processing, I tend to use [data.table package](https://rdatatable.gitlab.io/data.table/), which will be seen from now on.

```{r data-loading}
# Read files
old_quality_dt <-data.table()
for (f in old_files) {
    file_dt <- file.path(input_path, f) %>% fread()
    old_quality_dt <- rbind(old_quality_dt, file_dt)
}

new_quality_dt <- data.table()
for (f in new_files) {
    file_dt <- file.path(input_path, f) %>% fread()
    new_quality_dt <- rbind(new_quality_dt, file_dt)
}

str(old_quality_dt)
str(new_quality_dt)
```

# Preprocessament

Sent el format completament diferent, i considerant que volia totes les combinacions d'estacions, contaminants i temps, l'enfoc va ser:
1. Crear el dataset amb totes les combinacions desitjades, i separar-lo en dues parts (*vella* i *nova*).
2. Processar ambdós datasets separadament fins que es considerin estructurats, [tal com explica un gran gurú de R](https://vita.had.co.nz/papers/tidy-data.html).
3. Concatenar ambdós datasets per a crear el definitiu.

Així, el primer pas es va donar.

```{r quality_dt-structure}
stations_dt <- file.path(input_path, 'Qualitat_Aire_Estacions.csv') %>% fread()
stations_dt <- unique(stations_dt[, Codi_Contaminant:=NULL]) 
pollutants_dt <- file.path(input_path, 'Qualitat_Aire_Contaminants.csv') %>% fread()
pollutants_dt <- pollutants_dt[Desc_Contaminant %in% c('O3', 'NO2', 'PM10')]

date <- seq(limit_dates$min, limit_dates$max, 'days')
time <- 0:23
datetimes_dt <- CJ(date, time)

quality_dt <- merge(stations_dt[, key:=1], datetimes_dt[, key:=1], by='key', allow.cartesian=TRUE)
quality_dt <- merge(quality_dt[, key:=1], pollutants_dt[, key:=1], by='key', allow.cartesian=TRUE)
quality_dt[, key:=NULL]
names(quality_dt) <- names(quality_dt) %>% tolower()
quality_dt_old <- quality_dt[date < '2019-01-04']
quality_dt_new <- quality_dt[date >= '2019-01-04']
```

Després, les dades *velles* van ser processades.

```{r old-quality-preprocessing}
# In order, get day, convert hours to integer, and filter by pollutant
old_quality_dt[, date:=as.Date(gsub(pattern=' .*', replacement='', generat), format='%d/%m/%Y')]
old_quality_dt[, hora_o3:=as.integer(gsub(pattern='h', replacement='', hora_o3))]
old_quality_dt[, hora_no2:=as.integer(gsub(pattern='h', replacement='', hora_no2))]
old_quality_dt[, hora_pm10:=as.integer(gsub(pattern='h', replacement='', hora_pm10))]
old_quality_dt[, valor_o3:=as.numeric(gsub(pattern=' .*', replacement='', valor_o3))]
old_quality_dt[, valor_no2:=as.numeric(gsub(pattern=' .*', replacement='', valor_no2))]
old_quality_dt[, valor_pm10:=as.numeric(gsub(pattern=' .*', replacement='', valor_pm10))]

old_qual_o3_dt <- copy(old_quality_dt)
o3_code <- as.integer(pollutants_dt[Desc_Contaminant=='O3', .(Codi_Contaminant)])
old_qual_o3_dt <- old_qual_o3_dt[, .(codi_dtes, date, time=hora_o3, value=valor_o3, pollutant=o3_code)]
old_qual_o3_dt <- na.omit(old_qual_o3_dt)

old_qual_no2_dt <- copy(old_quality_dt)
no2_code <- as.integer(pollutants_dt[Desc_Contaminant=='NO2', .(Codi_Contaminant)])
old_qual_no2_dt <- old_qual_no2_dt[, .(codi_dtes, date, time=hora_no2, value=valor_no2, pollutant=no2_code)]
old_qual_no2_dt <- na.omit(old_qual_no2_dt)

old_qual_pm10_dt <- copy(old_quality_dt)
pm10_code <- as.integer(pollutants_dt[Desc_Contaminant=='PM10', .(Codi_Contaminant)])
old_qual_pm10_dt <- old_qual_pm10_dt[, .(codi_dtes, date, time=hora_pm10, value=valor_pm10, pollutant=pm10_code)]
old_qual_pm10_dt <- na.omit(old_qual_pm10_dt)

old_quality_dt <- rbind(old_qual_o3_dt, old_qual_no2_dt, old_qual_pm10_dt)
```

Les dades *noves* van ser processades també.

```{r new-quality-preprocessing}
# In order get date, from wide (hour) to values
names(new_quality_dt) <- names(new_quality_dt) %>% tolower()
new_quality_dt[, date:=as.Date(paste(as.character(any), as.character(mes), as.character(dia), sep='/'))]
value_cols <- names(new_quality_dt) %>% grep(pattern='^h[0-9][0-9]', value=TRUE)
new_quality_dt <- melt(new_quality_dt, id.vars=c('estacio', 'codi_contaminant', 'date'), measure.vars=value_cols, variable.name='time')
new_quality_dt[, time:=as.integer(gsub(pattern='h', replacement='', x=time)) -1]
```

Finalment, les dades van poder ser concatenades.

```{r dt-preprocessed}
# Join datasets
quality_dt_old <- merge(x=quality_dt_old, y=old_quality_dt, by.x=c('codi_dtes', 'date', 'time', 'codi_contaminant'),
		        by.y=c('codi_dtes', 'date', 'time', 'pollutant'), all.x=TRUE) 
quality_dt_new <- merge(x=quality_dt_new, y=new_quality_dt, by.x=c('estacio', 'date', 'time', 'codi_contaminant'),
		        by.y=c('estacio', 'date', 'time', 'codi_contaminant'), all.x=TRUE) 
quality_dt <- rbind(quality_dt_old, quality_dt_new)
setorder(quality_dt, estacio, codi_contaminant, date, time)
quality_dt[, datetime:=as.POSIXct(paste(date, time), format='%Y-%m-%d %H')]
quality_dt[, yearmonth:=format(date, '%Y%m')]
```

# Exploració

És mes fàcil entendre dades si aquestes estan agregades.
En primer lloc vaig fer una gràfica que mostrava valors mitjans per mes i contaminant, independentment de la estació on s'havia fet la mesura.

```{r monthly-eda-plot}
monthly_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, yearmonth)]
m <- ggplot(data=monthly_dt) + aes(x=yearmonth, y=value, color=desc_contaminant, group=desc_contaminant) + geom_line()
plot(m)
```

Es podia veure que no hi havia dades ni per febrer ni març de 2019, i que també els nivells de O3 són superiors a la resta de contaminants, i amb una variabilitat superior també.

Estava interessat amb l'evolució dels nivells de contaminació durant el dia, i per això també es va fer la gràfica de valors mitjans per hora i contaminant.

```{r hourly-eda-plot}
hourly_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, time)]
h <- ggplot(data=hourly_dt) + aes(x=time, y=value, color=desc_contaminant, group=desc_contaminant) + geom_line()
plot(h)
```

Es podia veure que hi havia una clara diferència entre contaminants: sent els nivells de PM10 i NO2 són relativament baixos després del migdia, els nivells de O3 són molt alts en comparació.
També hi havia la pregunta de si hi havia diferències entre estacions, independentment de l'hora.

```{r station-eda-plot}
stations_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, estacio)]
s <- ggplot(data=quality_dt) + aes(x=as.factor(estacio), y=value) + geom_boxplot() + facet_wrap(desc_contaminant~.)
plot(s)
```
En aquest cas hi havia valors extrems que no permetien veure les diferències. Es va comprovar que no hi havia cap patró, i es va procedir a excloure aquests valors per a la gràfica.

```{r station-eda-plot-no-outliers}
quality_dt[value > 1000]
s <- ggplot(data=quality_dt[value < 1000]) + aes(x=as.factor(estacio), y=value) + geom_boxplot() + facet_wrap(desc_contaminant~.)
plot(s)
```

No hi havia una clara diferència entre estacions, i per tant, entre districtes de Barcelona. El que va ser curiós és que els valors de PM10 estaven molt dispersats per sobre del percentil 75.
Havent vist estacions sense pràcticament dades, vaig decidir mantenir només les combinacions de estacions i contaminant que contenen valors.

```{r exclude-combinations-without-values}
perc_na_dt <- quality_dt[, .(perc_nas=sum(is.na(value))/.N), by=.(desc_contaminant, estacio)]
include_comb_dt <- perc_na_dt[perc_nas < 0.8, .(desc_contaminant, estacio)]  
quality_dt <- merge(quality_dt, include_comb_dt, by=c('desc_contaminant', 'estacio'))
```

Considerant valors extrems, i tenint en compte que els mètodes de sèries temporals en general no gestionen valors no definits (NA), es va decidir intentar corregir ambdós problemes.

```{r remove-outliers-nas}
# Replace NA values and outliers - modelling purposes
quality_dt[, original_value:=value]
quality_dt[, value:=as.numeric(tsclean(as.ts(original_value))), by=.(desc_contaminant, estacio)]
quality_dt[, imputed_values:=0]
quality_dt[is.na(original_value) | original_value != value, imputed_values:=1]
```

Després d'una búsqueda ràpida, vaig trobar [estàndards definits per la Organització Mundial de la Salut](https://apps.who.int/iris/bitstream/handle/10665/69477/WHO_SDE_PHE_OEH_06.02_eng.pdf) 
relatiu als 3 contaminants considerats. Així, podria ser molt útil visualitzar els nivells actuals en comparació amb els estàndars.

```{r standards-comparison}
# Compare to standards
# PM10 standard: 50/24hour
# O3 standards: 100/8hour
# NO2 standards: 200/hour

standards_pm10 <- 50
standards_o3 <- 100
standards_no2 <- 200

daily_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, date)]
quality_dt[, hour_range:=NaN]
quality_dt[between(time, 0, 7), hour_range:=1]
quality_dt[between(time, 8, 15), hour_range:=2]
quality_dt[between(time, 16, 23), hour_range:=3]
daily_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, date)]
hrange_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, date, hour_range)]

pm10_p <- ggplot(data=daily_dt[desc_contaminant=='PM10']) + aes(x=date, y=value) + geom_line() +
	geom_hline(yintercept=standards_pm10, color='red')
o3_p <- ggplot(data=hrange_dt[desc_contaminant=='O3']) + aes(x=date, y=value) + geom_line() +
	geom_hline(yintercept=standards_o3, color='red')
no2_p <- ggplot(data=quality_dt[desc_contaminant=='NO2']) + aes(x=date, y=value) + geom_line() +
	geom_hline(yintercept=standards_no2, color='red')
plot(pm10_p)
plot(o3_p)
plot(no2_p)
```

# Modelització

Predir els nivells de contaminació pot proveïr d'informació als ajuntaments, amb la qual cosa podria mitigar nivells alts de predicció, per exemple.

En general, predicció de sèries temporals pot no ser molt acurada per prediccions a llarg termini; per aquest motiu només els valors de octubre de 2019 són els que es prediran. 

Per simplificació i motius computacionals, dades diàries (mitjana de valors per dia) enlloc de dades horàries va ser utilitzat per a l'anàlisi.

Per tal de saber el potencial comportament del model sobre noves dades, va ser necessari de fer una prova a passat prèviament.
Això implica que podíem entrenar un model de juny del 2018 a novembre del 2018, predint valors de desembre de 2018, i així poder comparar-les amb la realitat.
El model que millor funcioni a passat seria l'utilitzat a futur.

En aquest cas, els dos models a comparar van ser: [ARIMA](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) i [Exponential Smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing).
La mètrica per a comparar varis models era [SMAPE](https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error), la qual considera diferències relatives.


```{r model-comparison}
gen_model_dates <- list(train_ini=as.Date('2018/06/01'), train_end=as.Date('2018/11/30'), test_end=as.Date('2018/12/31'))
daily_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, date, estacio)]

test_models <- function(dt, model_dates) {
	train_ts <- dt[between(date, model_dates$train_ini, model_dates$train_end), value]
	test_ts <- dt[between(date, model_dates$train_end + 1, model_dates$test_end), value]
	fc_periods <- as.numeric((model_dates$test_end) - (model_dates$train_end + 1) + 1)
	pred_autoarima <- auto.arima(train_ts) %>% forecast(fc_periods)
	pred_ets <- ets(train_ts) %>% forecast(fc_periods)
	results_dt <- data.table(actuals=test_ts, pred_autoarima=as.numeric(pred_autoarima$mean), 
						  pred_ets=as.numeric(pred_ets$mean))
	return(results_dt)
}

test_dt <- daily_dt[, test_models(dt=.SD, model_dates=gen_model_dates), by=.(desc_contaminant, estacio)]
smape_autoarima <- smape(test_dt[, actuals], test_dt[, pred_autoarima])
smape_ets <- smape(test_dt[, actuals], test_dt[, pred_ets])
cat('Smape Autoarima', smape_autoarima)
cat('Smape ETS', smape_ets)
```

Tenint en compte que el model d'Autoarima va ser el de menor mètrica de les dues opcions, va ser el model escollit per a ser utilitzat amb les dades més recents.

```{r predictions}
apply_model_dates <- list(train_ini=as.Date('2019/04/01'), train_end=as.Date('2019/09/30'), test_end=as.Date('2019/10/31'))

apply_models <- function(dt, model_dates) {
	train_ts <- dt[between(date, model_dates$train_ini, model_dates$train_end), value]
	fc_periods <- as.numeric((model_dates$test_end) - (model_dates$train_end + 1) + 1)
	pred_autoarima <- auto.arima(train_ts) %>% forecast(fc_periods)
	pred_dates <- seq(model_dates$train_end + 1, model_dates$test_end, 'days')
	pred_dt <- data.table(value=as.numeric(pred_autoarima$mean), date=pred_dates)
	return(pred_dt)
}

predictions_dt <- daily_dt[, apply_models(dt=.SD, model_dates=apply_model_dates), by=.(desc_contaminant, estacio)]
daily_dt <- rbind(daily_dt, predictions_dt)
```

# Aplicació

L'últim pas va ser preparar el dataset que va ser utilitzat per a la aplicació interactiva. Variables importants eren les següents:
- Contaminant
- Districte
- Districte i coordenades de l'estació de contaminació
- Data
- Estàndars de la OMS
- Mètriques

Per motius d'experiència d'usuari, vaig pensar que seria millor visualitzar dades a nivell mensual. La pregunta llavors era, quines mètriques s'haurien de considerar per a aquest nou dataset.

En aquest case, les següents tres van ser les mètriques definides:
- Mitjana: Valor mitjà, tenint en compte valors mitjans per dia (mitjana de mitjanes).
- Màxima: Valor màxim, tenint en compte valors mitjans per dia (màxim de mitjanes).
- Percentatge per sobre dels estàndars: Percentatge de dies el màxim dels quals està per sobre dels estàndars de la OMS.

**És important considerar que les comparacions amb els estàndards mencionats poden no ser completament justes.
Això rau en què els estàndars de la OMS consideren un rang temporal diferent per a cada contaminant, mentre que les mètriques són calculades com a mitjanes, màxims i percentatges de mitjanes diàries.**

```{r final-processing}
# Data for app - select and save dataset
stations_dt <- unique(quality_dt[, .(estacio, nom_districte, latitud, longitud)])
daily_dt <- merge(daily_dt, stations_dt, by=c('estacio'))
daily_dt <- daily_dt[, estacio:=NULL]
old_names <- c('desc_contaminant', 'nom_districte', 'latitud', 'longitud')
new_names <- c('pollutant', 'district', 'latitude', 'longitude')
setnames(daily_dt, old_names, new_names)

units <- unique(pollutants_dt[, .(units=Unitats)])
standards_dt <- data.table(pollutant=c('PM10', 'NO2', 'O3'), standards=c(standards_pm10, standards_no2, standards_o3))
standards_dt[, units:=units]
daily_dt <- merge(daily_dt, standards_dt, by='pollutant')
daily_dt[, yearmonth:=format(date, '%Y%m')]

monthly_dt <- daily_dt[, .(mean_value=round(mean(value), 2), max_value=round(max(value), 2), perc_above_st=sum(value>standards)/.N), 
			by=.(pollutant, yearmonth, district, latitude, longitude, standards, units)] 

fwrite(monthly_dt, file.path(output_path, 'app_data.csv'))
```

[Shiny](https://shiny.rstudio.com/) és un paquet de R per a desenvolupar informes interactius a partir de codi R. Aquest va ser el paquet utilitzat per a la aplicació final.
L'estructura genèrica d'una aplicació shiny és la següent:
- global.R: Carrega les dades i les processa si és necessari.
- server.R: Genera objectes d'interacció amb les dades.
- ui.R: Crea la interfície de l'usuari.

El codi global.R és el que val la pena detallar.

En la aplicació, hi ha cercles que representen estacions, i ambdós radis i colors són parametritzables; per aquest motiu dues columnes auxiliars van ser generades.

La lògica de colors es va definir de la manera següent:
- Si la mètrica *Max* és superior als estàndards, l'estació apareixerà de color vermell.
- Si la mètrica *Max* és superior o igual al 80% dels estàndards, l'estació apareixerà de color groc.
- Si la mètrica *Max* és inferior al 80% dels estàndars, l'estació apareixerà de color verd.

L'intuició dels radis és que els més grans indicaven valors més extrems (fossin aquests bons o dolents). La lògica concreta era la següent:
- Color verd:
	- El radi màxim es pot veure quan els valors van del 0 al 50% dels estàndars.
	- El radi decreixerà linealment quan els valors s'apropin al 80% dels estàndars, convergint al radi mínim.
- Color groc: El radi és constant i el seu valor és el mínim.
- Color vermell:
	- El radi creixerà linealment de mínim a màxim quan els valors s'apropin als 200% dels estàndars.
	- El radi màxim es podrà veure quan els valors siguin 200% dels estàndars o superiors.

```{r global-app }
data_path <- 'data/output'
dt <- file.path(data_path, 'app_data.csv') %>% fread()

# Add features to data for map usage
colour_ratio <- 0.8

dt[, risk_colour:='green']
dt[max_value >= colour_ratio * standards, risk_colour:='yellow']
dt[max_value > standards, risk_colour:='red']

min_radius <- 20
max_radius <- 40
min_ratio <- 0.5
max_ratio <- 2

dt[, risk_radius:=min_radius]
dt[risk_colour=='red', risk_radius:=min_radius * (max_value / standards)]
dt[risk_colour=='red' & max_value > max_ratio * standards, risk_radius:=max_radius]
dt[risk_colour=='green', risk_radius:=min_radius * (standards / max_value)]
dt[risk_colour=='green' & max_value < min_ratio * standards, risk_radius:=max_radius]

```

La resta del codi pel desenvolupament de l'aplicació no està detallat en l'anàlisi. Com ja s'ha comentat, es pot trobar en el repositori especificat.

# Finalització

El principal objectiu d'aquest anàlisi ha estat mostrar el potencial de:
- Tecnologies de codi obert
- Dades obertes
- Modelització de dades
- Eines de visualització

Usant-les apropiadament, és possible tenir un impacte positiu en la qualitat de vida de les persones.

R dataviz2019_eng.Rmd => dataviz_eng.Rmd +137 -27
@@ 1,43 1,53 @@
---
title: 'Analyzing pollutant levels from Barcelona'
title: 'Analysis of pollution levels in Barcelona'
output:
  html_document
---

# Introduction

This analysis has been made in order to participate in the [World Dataviz Challenge 2019](https://opendata-ajuntament.barcelona.cat/ca/data-viz-kobe-home).
There are huge environmental concerns nowadays, so I decided that it could be useful to try to get insights, predict and visualize pollution levels in the city of Barcelona.
[R](https://www.r-project.org/), one of the main programming languages used in data science (and my personal favourite), has been the chosen one for this development.
This analysis was made in order to participate in the [World Dataviz Challenge 2019](https://opendata-ajuntament.barcelona.cat/en/data-viz-kobe-home).
Given the growing environmental concerns nowadays, and the explicit recommendation of the challenge to address pollution as a topic, the choice was pretty clear.
[R](https://www.r-project.org/), one of the main programming languages used in data science (and my personal favourite), was the chosen one for this development.

This document will cover, very briefly and not thoroughly, the main steps that can be found in a data science pipeline:
1. Loading: Load the data
The analysis consisted in the following steps:
1. Loading: Load pollution data.
2. Preprocessing: Preprocess the data until it is structured enough for analyze.
3. Exploration data analysis: Get insights from the data, usually using plots.
4. Modelling & predicting: Train and apply a model to predict future outcome.
5. Application: Develop an application that embeds all previous work.
3. Exploration: Get insights from pollution trends and characteristics using plots.
4. Modelling: Train and apply a model to predict future pollution levels.
5. Application: Develop an application that show pollution levels evolution by station.

# Data & packages loading
The source code for the analysis, documents and application can be found at [sourcehut](https://git.sr.ht/~thrill_seeker/dataviz-challenge-2019).
The application can be seen at [shiny portal](x).

First, of course it was required to load the data. They were retrieved from the following links:
- [Air data quality](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/qualitat-aire-detall-bcn)
- [Pollution measure stations](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/qualitat-aire-estacions-bcn)
- [Pollutants](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/contaminants-estacions-mesura-qualitat-aire )
# Loading

First, it was required to load the data. They were retrieved from the following links:
- [Air data quality](https://opendata-ajuntament.barcelona.cat/data/en/dataset/qualitat-aire-detall-bcn)
- [Pollution measure stations](https://opendata-ajuntament.barcelona.cat/data/en/dataset/qualitat-aire-estacions-bcn)
- [Pollutants](https://opendata-ajuntament.barcelona.cat/data/en/dataset/contaminants-estacions-mesura-qualitat-aire )

Althought there is the possibility to use an API directly to retrieve the data, in this case I decided to just download the .csv files.

The required packages for development are loaded.
The required packages for development were loaded.

```{r package-loading}
library(data.table)
library(dplyr)
library(forecast)
library(ggplot2)
library(leaflet)
library(Metrics)
library(scales)
library(shiny)
library(shinydashboard)
library(shinyWidgets)
```
As specified in the Air Data Quality URL, dataset structure have changed completely from April 2019 onwards.
I decided then to define two different loading parts, the one that gets the data from the June 2018 (first file available) until March 2019 (the last with the same structure),
and the one that gets data from April 2019 to September 2019 (the last file at that moment). Auxiliary variables were also defined.

As specified in the Air Data Quality URL, dataset structure had changed completely since April 2019.
I decided then to define two different loading parts:
1. The *old* one that gets the data from the June 2018 (first file available) until March 2019 (the last with the same structure).
2. The *new* one that gets data from April 2019 to September 2019 (the last file at that moment). Auxiliary variables were also defined.

```{r variable-loading}
# Define variables


@@ 74,11 84,11 @@ str(old_quality_dt)
str(new_quality_dt)
```

# Data preprocessing
# Preprocessing

Being the format completely different, and given that I wanted to have all combinations of station, pollutant, and times, the approach followed was:
1. Create a dataset will all the desired combinations, and separate it into two parts (old and new structure).
2. Process old data and new data separatly for both to be considered tidy, [as explained by a great R guru](https://vita.had.co.nz/papers/tidy-data.html).
1. Create a dataset will all the desired combinations, and separate it into two parts (*old* and *new*).
2. Process both datasets separately for them to be considered tidy, [as explained by a great R guru](https://vita.had.co.nz/papers/tidy-data.html).
3. Concatenate both datasets to create the definite one.

Then, the first step was taken.


@@ 156,7 166,7 @@ quality_dt[, datetime:=as.POSIXct(paste(date, time), format='%Y-%m-%d %H')]
quality_dt[, yearmonth:=format(date, '%Y%m')]
```

# Exploratory data analysis
# Exploration

It is easier to understand data if they are aggregated.
I plotted mean values by month and pollutant, regardless of the station in which they are measured.


@@ 203,7 213,7 @@ include_comb_dt <- perc_na_dt[perc_nas < 0.8, .(desc_contaminant, estacio)]
quality_dt <- merge(quality_dt, include_comb_dt, by=c('desc_contaminant', 'estacio'))
```

Considering outliers, and taking into account that forecasting methods in general do not handle NA values, it was decided to correct both problems.
Considering outliers, and taking into account that time series forecasting methods in general do not handle NA values, it was decided to correct both problems.
In this case, it was done using a function called tsclean, from the forecast package.

```{r remove-outliers-nas}


@@ 215,7 225,7 @@ quality_dt[is.na(original_value) | original_value != value, imputed_values:=1]
```

After a quick search, I found [standards defined by the World Health Organization](https://apps.who.int/iris/bitstream/handle/10665/69477/WHO_SDE_PHE_OEH_06.02_eng.pdf) 
related to all 3 pollutants considered: O3, NO2 and PM10. Then, it could be very useful to visualize their current levels with the standards defined.
related to all 3 pollutants considered: O3, NO2 and PM10. Then, it could be very useful to visualize their current levels compared to the standards defined.

```{r standards-comparison}
# Compare to standards


@@ 248,7 258,7 @@ plot(no2_p)

# Modelling

Predict future pollution levels can provide information to city councils, which can be used to mitigate high-predicted levels, for example.
Predict future pollution levels can provide information to city councils, which can be used to mitigate high-predicted pollution levels, for example.

In general, time series forecasting may not be really accurate for long-term predictions; this is why just October 2019 values will be predicted.
Taking into account that there were no data for February and March 2019, the idea was to use from April 2019 to September 2019 to train the model.


@@ 284,7 294,7 @@ cat('Smape Autoarima', smape_autoarima)
cat('Smape ETS', smape_ets)
```

As Autoarima had the lowest smape of the two options, was decided to be applied to recent data.
As Autoarima had the lowest smape of the two options, it was the chosen model to be applied to recent data.

```{r predictions}
apply_model_dates <- list(train_ini=as.Date('2019/04/01'), train_end=as.Date('2019/09/30'), test_end=as.Date('2019/10/31'))


@@ 302,4 312,104 @@ predictions_dt <- daily_dt[, apply_models(dt=.SD, model_dates=apply_model_dates)
daily_dt <- rbind(daily_dt, predictions_dt)
```

# Prepare application
# Application

The last step was to prepare the dataset that would be used for the interactive application. Important features where the following:
- Pollutant
- District
- Pollution station district and coordinates
- Date
- WHO Standards
- Metrics

For user experience purposes, I thought it would be better to visualize monthly data. The question then was, which metrics should be taken into account for this new dataset.

In this case, three were the metrics chosen:
- Mean: Mean value, taking into account mean values per day (mean of means).
- Max: Maximum value, taking into account mean values per day (max of means).
- Percentage above standards: Percentage of days whose maximum value is above the defined WHO standards.

**It is important to note that comparisons to WHO standards may not be completely fair. 
This is because WHO standards consider different time ranges for mean values, whereas this metrics are computed as monthly means, maximums and percentages from daily means.**

```{r final-processing}
# Data for app - select and save dataset
stations_dt <- unique(quality_dt[, .(estacio, nom_districte, latitud, longitud)])
daily_dt <- merge(daily_dt, stations_dt, by=c('estacio'))
daily_dt <- daily_dt[, estacio:=NULL]
old_names <- c('desc_contaminant', 'nom_districte', 'latitud', 'longitud')
new_names <- c('pollutant', 'district', 'latitude', 'longitude')
setnames(daily_dt, old_names, new_names)

units <- unique(pollutants_dt[, .(units=Unitats)])
standards_dt <- data.table(pollutant=c('PM10', 'NO2', 'O3'), standards=c(standards_pm10, standards_no2, standards_o3))
standards_dt[, units:=units]
daily_dt <- merge(daily_dt, standards_dt, by='pollutant')
daily_dt[, yearmonth:=format(date, '%Y%m')]

monthly_dt <- daily_dt[, .(mean_value=round(mean(value), 2), max_value=round(max(value), 2), perc_above_st=sum(value>standards)/.N), 
			by=.(pollutant, yearmonth, district, latitude, longitude, standards, units)] 

fwrite(monthly_dt, file.path(output_path, 'app_data.csv'))
```

**It is important to note that comparisons to WHO standards may not be completely fair. 
This is because WHO standards consider different time ranges for mean values, whereas this metrics are computed as monthly means, maximums and percentages from daily means.**

[Shiny](https://shiny.rstudio.com/) is an R package to develop interactive dashboards from R code. This was the package used for the final application.
The structure of a generic shiny application is the following:
- global.R: Loads the data and process them if necessary.
- server.R: Generates build objects that interact with the data.
- ui.R: Generates the user interface.

The code global.R is the one worth detailing.

In the application, there are circles representing stations, and both radius and colour are customizable, that is why two auxiliar columns will be generated.

Colour logic will be defined as follows:
- If *Max* metric is greater than standards, station colour will be red.
- If *Max* metric is greater or equal than 80% of standards, station colour will be yellow.
- If *Max* metric is lower than 80% of standards, station colour will be green.

Radius intuition is that bigger radius means more extreme values (either good or bad). The specific logic is the following:
- Green colour:
	- Maximum radius can be seen when values range from 0 to 50% of standards.
	- Radius will linearly decrease as values approaches 80% of standards, converging then to minimum radius.
- Yellow colour: Radius is constant and its value set to minimum.
- Red colour:
	- Radius will linearly increase from minimum to maximum as values approach 200% of standards.
	- Maximum radius can be seen when values are from 200% of standards onwards.

```{r global-app }
data_path <- 'data/output'
dt <- file.path(data_path, 'app_data.csv') %>% fread()

# Add features to data for map usage
colour_ratio <- 0.8

dt[, risk_colour:='green']
dt[max_value >= colour_ratio * standards, risk_colour:='yellow']
dt[max_value > standards, risk_colour:='red']

min_radius <- 20
max_radius <- 40
min_ratio <- 0.5
max_ratio <- 2

dt[, risk_radius:=min_radius]
dt[risk_colour=='red', risk_radius:=min_radius * (max_value / standards)]
dt[risk_colour=='red' & max_value > max_ratio * standards, risk_radius:=max_radius]
dt[risk_colour=='green', risk_radius:=min_radius * (standards / max_value)]
dt[risk_colour=='green' & max_value < min_ratio * standards, risk_radius:=max_radius]
```
The rest of code for shiny deployment is not detailed in this analysis. As has been already said, it can be found in the specified repository.

# Wrap up

The main objective of the analysis is to show the potential of using:
- Open source technologies
- Open data
- Data modelling
- Visualization tools

Using them appropriately, it is possible to have an impact on people's quality of life.