This visualization project is based on real open data. You can access all available datasets at: Câmara Municipal de Lisboa

library(tidyverse)
library(leaflet)
library(mapview)
library(htmlwidgets)
library(viridis) 
library(scales)  
library(plotly)
library(kableExtra)
library(shiny)
library(leaflet.minicharts)
library(RColorBrewer)
library(magick)

data1 <- read_csv("QA_2023_Q1.csv")
data2 <- read_csv("QA_2023_Q2.csv")
data3 <- read_csv("QA_2023_Q3.csv")
data4 <- read_csv("QA_2021_Q2.csv")
data5 <- read_csv("QA_2021_Q3.csv")
data6 <- read_csv("QA_2022_Q1.csv")
data7 <- read_csv("QA_2022_Q2.csv")
data8 <- read_csv("QA_2022_Q3.csv")
data9 <- read_csv("QA_2024_Q1.csv")

data_list <- mget(paste0("data", 1:9))

QAall <- bind_rows(data_list) |>
  filter(VALOR <1000 & VALOR >0 & PARAMETRO != "NO") |>    # exclude outliers and negative values; and NO (unimportant pollutant)
  dplyr::select(c(DTM_UTC, PARAMETRO, LOCAL, LATITUDE, LONGITUDE, COR_NIVEL, VALOR)) %>%    # exclude unimportant columns
  dplyr::rename(
    moment = DTM_UTC,  #datetime column
    station = LOCAL   #name of measuring stations
  ) |> 
  mutate(
    date = as.Date(moment),
    day = day(moment),
    year = year(moment),
    month = month(moment),
    hour = hour(moment),
    weekday =wday(moment, label = TRUE, locale = "en_US")
  ) %>% 
  arrange(date)

Pollutants analyzed

Figure 1

There are 79 air quality monitoring stations in Lisbon, although not all stations monitor all pollutants. This interactive map shows the locations of the stations and the pollutants analyzed.

#Table with stations and their localization and pollutants analyzed
stationvalues <- QAall %>%
  pivot_wider(names_from = PARAMETRO, values_from = VALOR) %>%
  group_by(station) %>%
  summarise(
    long = first(LONGITUDE),
    lat = first(LATITUDE),
    CO = ifelse(mean(CO, na.rm = TRUE) == 0 | is.na(mean(CO, na.rm = TRUE)),0,1),
    NO2 = ifelse(mean(NO2, na.rm = TRUE)  == 0 | is.na(mean(NO2, na.rm = TRUE)),0,1),
    O3 = ifelse(mean(O3, na.rm = TRUE)  == 0 | is.na(mean(O3, na.rm = TRUE)),0,1),
    PM25 = ifelse(mean(PM25, na.rm = TRUE) == 0 | is.na(mean(PM25, na.rm = TRUE)),0,1),
    PM10 = ifelse(mean(PM10, na.rm = TRUE) == 0 | is.na(mean(PM10, na.rm = TRUE)),0,1),
    SO2 = ifelse(mean(SO2, na.rm = TRUE) == 0 | is.na(mean(SO2, na.rm = TRUE)),0,1),
  ) %>%
  mutate(station_code = row_number())

pollutant_colors <- c("#FF0000", "purple", "yellow", "#00FF00", "#0000FF", "#00FFFF")


map <- leaflet() %>%
  addTiles() %>%
  addMinicharts(
    stationvalues$long, stationvalues$lat,
    type = "pie",
    chartdata = dplyr::select(stationvalues, NO2, O3, CO, PM10, PM25, SO2),
    colorPalette = pollutant_colors,
    popup = popupArgs(
      showTitle = F,
      showValues = F)
  ) |> 
  addCircleMarkers(
    lng = stationvalues$long, 
    lat = stationvalues$lat, 
    label = stationvalues$station,
    color = "#00000000"      
  )
map

Figure 1. Air Quality Monitoring: 79 Stations Across Lisbon Measuring Key Pollutants.

Figure 2

This bar plot shows the annual NO2 mean levels in 2023 for each station. Colors (green for < 40 µg/m³ and yellow to red for higher levels) are based on the 2005 guidelines If the new 2021 guidelines were used, all stations would have exceeded the limit (10 µg/m³).

stationvalues23 <- QAall %>%
  filter(year == 2023) |> 
  pivot_wider(names_from = PARAMETRO, values_from = VALOR) %>%
  group_by(station) %>%
  summarise(
    long = first(LONGITUDE),
    lat = first(LATITUDE),
    NO2 = mean(NO2, na.rm = TRUE),
    O3 = mean(O3, na.rm = TRUE),
    PM25 = mean(PM25, na.rm = TRUE),
    PM10 = mean(PM10, na.rm = TRUE),
    SO2 = mean(SO2, na.rm = TRUE)
  )


# Function to assign colors based on NO2 values
assign_colors <- function(value) {
  if (value <= 40) {
    return(scales::rescale(value, to = c(0, 1)) * 120)  # Scale for green to yellow
  } else {
    return(scales::rescale(value, to = c(0, 1)) * 120 + 120)  # Scale for yellow to red
  }
}

# Assign colors using a color palette for values <= 40 and > 40
stationvalues23 <- stationvalues23 |>
  filter(!is.nan(NO2)) |>
  arrange(NO2) |>
  mutate(
    station = factor(station, levels = unique(station)),
    color = case_when(
      NO2 <= 40 ~ scales::col_numeric(palette = c("green", "yellow"), domain = NULL)(NO2),
      NO2 > 40  ~ scales::col_numeric(palette = c("yellow", "red"), domain = NULL)(NO2)
    )
  )

p <- ggplot(stationvalues23) +
  geom_col(aes(x = NO2, y = fct_infreq(station), fill = color)) +
  geom_vline(xintercept = 40, color = "red", linetype = "dashed", size = 1) +
  geom_vline(xintercept = 10, color = "black", linetype = "dashed", size = 1) +
  scale_fill_identity() +  # Use the color column as fill color
  labs(title =  "Anual mean NO2 values per station", x = "NO2 (ug/m3)") +
  scale_x_continuous(expand = c(0, 0)) +
  theme(
    legend.position = "none",
    panel.background = element_rect(fill = "white"),
    panel.border = element_blank(),
    panel.grid.major = element_line(linewidth = 0.5, linetype = 'solid', colour = "gray"),
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11),
    axis.title = element_text(size = 12),
    axis.title.y = element_blank(),
    axis.title.x = element_text(margin = margin(t = 10))
  ) 
p
Figure 2. Annual mean NO2 levels across stations in 2023. Note the older 2005 limit (40 ug/m³, dashed red) vs. the newer 2021 limit (10 ug/m³, dashed black) from the WHO air quality guidelines.

Figure 2. Annual mean NO2 levels across stations in 2023. Note the older 2005 limit (40 ug/m³, dashed red) vs. the newer 2021 limit (10 ug/m³, dashed black) from the WHO air quality guidelines.

Figure 3

This heatmap shows the concentration of pollutants by weekday and hour of the day, highlighting the correlation with rush hours (morning and afternoon). Despite an overall improvement in pollutant concentrations since 2021, levels remain high.

save_heatmap <- function(data, parameter, filename) {
  heat <- data |>
    filter(year <= 2023, PARAMETRO == parameter) |>
    mutate(weekday = factor(weekday, levels = c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"))) |>
    group_by(PARAMETRO, year, weekday, hour) |>
    summarise(meanvalue = mean(VALOR, na.rm = TRUE), .groups = 'keep')
  
  p <- heat |>
    ggplot(aes(x = weekday, y = hour, fill = meanvalue)) + 
    geom_tile(color = "white", size = 0.1) +
    facet_wrap(~ year, ncol = 3) + 
    scale_fill_stepsn(name = "ug/m3", n.breaks = 9, colours = viridis::inferno(9)) +
    scale_y_continuous(expand = c(0, 0), breaks = c(0, 6, 12, 18, 23)) +
    coord_equal() + # make tiles squares
    labs(x = NULL, y = NULL, title = paste("Mean level of", parameter, "per weekday and time of day")) +
    theme(plot.title = element_text(size = 16, hjust = 0),
          axis.ticks = element_blank(),
          axis.text.y = element_text(size = 11),
          axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 0.9),
          legend.title = element_text(size = 13),
          legend.text = element_text(size = 11),
          strip.text = element_text(size = 13)
    )
  
  ggsave(filename, plot = p, width = 8, height = 5, dpi = 100)
}

parameters <- unique(QAall$PARAMETRO)
filenames <- paste0(parameters, ".png")
for (i in seq_along(parameters)) {
  save_heatmap(QAall, parameters[i], filenames[i])
}

# Create a GIF from the saved images
images <- image_read(filenames)
animation <- image_animate(images, fps = 0.25) 
animation
Figure 3. Heatmap of pollutant concentrations (by hour and weekday).
Figure 3. Heatmap of pollutant concentrations (by hour and weekday).

Figure 4

This interactive timeline shows the daily mean levels of pollutants since 2021. Some days have missing data (represented by a straight line between dates, especially in April 2023).

summarised_data <- QAall %>%
  filter(PARAMETRO != "CO") |>  # CO has other units (mg/m3), and all measurements were low.
  group_by(date, PARAMETRO) %>%
  summarise(VALOR = mean(VALOR), .groups = 'drop')

# Custom date label function
custom_date_labels <- function(x) {
  ifelse(format(x, "%m") == "01" | (format(x, "%m") == "06" & format(x, "%Y") == "2021"), 
         format(x, "%b-%Y"), 
         format(x, "%b"))
}

timeline <- ggplot(data = summarised_data, aes(x = date, y = VALOR, color = PARAMETRO)) +
  geom_line() +
  scale_color_manual(values = c("NO2" = "blue", "PM10" = "red", "PM25" = "green", "O3" = "brown", "SO2" = "orange")) +
  geom_hline(yintercept = 15, colour = "green", linetype = "dashed") +
  geom_hline(yintercept = 40, colour = "orange", linetype = "dashed") +
  geom_hline(yintercept = 45, colour = "red", linetype = "dashed") +
  labs(title = "Daily mean level of pollutants", y = "Concentration (ug/m3)", color= "Pollutant")+
  geom_hline(yintercept = 25, colour = "blue", linetype = "dashed") +
  scale_x_date(date_breaks = "1 month", labels = custom_date_labels) +
  theme(
    panel.background = element_rect(fill = "white"),
    panel.border = element_blank(),
    panel.grid.major = element_line(linewidth = 0.5, linetype = 'solid', colour = "gray"),
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.text.x = element_text(size = 11, angle = 45, hjust = 0.5, vjust = 0.5),
    axis.text.y = element_text(size = 11),
    axis.title = element_text(size = 12),
    axis.title.x = element_blank()
    )

ggplotly(timeline,
         width = 1200,
         modebar = list(displayModeBar = FALSE)
        )

Figure 4. Mean concentration of pollutants for each day (August 2021 - January 2024). Dashed lines represent the 2021 WHO air quality guidelines daily limits.

Figure 5

This area chart shows the proportion of pollutants for each date from 2021 to 2023. NO2 and O3 account for almost 75% of the contaminants.

stacked <- QAall |> 
  mutate(date = as.Date(paste(year, month, "01", sep = "-"))) |> 
  group_by(date, PARAMETRO) |>
  summarise(value = mean(VALOR, na.rm = TRUE))



custom_date_format <- function(x) {
  format_date <- function(date) {
    ifelse(format(date, "%m") == "01",
           format(date, "%b-%y"), 
           format(date, "%b"))
  }
  months <- format_date(x)
  months_capitalized <- tools::toTitleCase(months)
  return(months_capitalized)
}

ggplot(stacked, aes(x = date, y = value, fill = PARAMETRO)) +
  geom_area(position = "fill", stat = "identity", colour="black") +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Proportion of Pollutants by date (2021-2024)", y = "Concentration ug/m3") +
  scale_x_date(labels = custom_date_format, date_breaks = "1 month", expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0))+
  theme(legend.title = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title.x = element_blank())
Figure 5. Area chart showing the proportion of pollutants by date.

Figure 5. Area chart showing the proportion of pollutants by date.

Figure 6

This map shows the main pollutant for each station. Since each pollutant has different concentration ranges, the data was normalized.

maxpol_normalized <- QAall %>%
  filter(year == 2023 & PARAMETRO != "CO") |> 
  group_by(PARAMETRO) %>%
  mutate(normalized_valor = (VALOR - mean(VALOR, na.rm = TRUE)) / sd(VALOR, na.rm = TRUE),
         long = first(LONGITUDE),
         lat = first(LATITUDE)) %>%
  ungroup()

# Calculate both raw and normalized means and identify the predominant pollutant
result <- maxpol_normalized %>%
  group_by(station, PARAMETRO) %>%
  summarise(
    raw_mean_valor = mean(VALOR, na.rm = TRUE),                   
    normalized_mean_valor = mean(normalized_valor, na.rm = TRUE),
    long = first(LONGITUDE),
    lat = first(LATITUDE),
  ) %>%
  arrange(station, desc(raw_mean_valor)) %>%                       
  mutate(
    predominant_raw = ifelse(row_number() == 1, PARAMETRO, NA),   
    predominant_normalized = ifelse(normalized_mean_valor == max(normalized_mean_valor), PARAMETRO, NA)
  ) %>%                                                          
  filter(!is.na(predominant_raw) | !is.na(predominant_normalized)) %>% 
  ungroup() %>%
  select(station, predominant_raw, raw_mean_valor, predominant_normalized, normalized_mean_valor, lat, long)

palette <- colorFactor(
  palette = "Set1",   
  domain = result$predominant_normalized
)

result <- result |> 
  filter(!is.na(predominant_normalized))

# Map with different colors for each main pollutant
lisbon_map <- leaflet(result) %>%
  addTiles() %>%
  setView(lng = -9.1393, lat = 38.7223, zoom = 12) %>%
  addCircleMarkers(
    ~long, 
    ~lat, 
    color = ~palette(predominant_normalized),  
    popup = ~predominant_normalized, 
    label = ~station,
    radius = 6,  
    fillOpacity = 0.8,  
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright", 
    pal = palette, 
    values = ~predominant_normalized,
    title = "Main Pollutant",
    opacity = 1
  )

lisbon_map

Figure 6. Main pollutant per station (normalized).

Figure 7

This plot represents the mean levels for each month of the year. Note that O₃ and SO₂ levels are higher in the warmer months, while PM₁₀ and PM₂.₅ levels are higher during winter. NO₂ levels decrease sharply in August due to holidays and reduced traffic in the city.

monthly <- QAall |> 
  filter(PARAMETRO != "CO") |> 
  group_by(month, PARAMETRO) |> 
  summarise(value = mean(VALOR, na.rm = TRUE)) |> 
  mutate(month = factor(month, 
                        levels = 1:12, 
                        labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", 
                                   "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")))

ggplot(monthly, aes(x = month, y = value, colour = fct_reorder2(PARAMETRO, month, value), group = PARAMETRO)) +
  geom_line(size = 1.2) + 
  labs(x = "Month", y = "Mean concentration (ug/m3)", title = "Mean Values by Pollutant and month", 
       colour = "Pollutant") +
  scale_x_discrete(expand = c(0,0)) +
  theme_bw()
Figure 7. Mean Values by month.

Figure 7. Mean Values by month.

Figure 8

This line plot represents the daily CO levels over time. Although CO concentration has risen recently, the daily limit according to WHO 2021 guidelines is 4 mg/m³.

custom_date_labels2 <- function(x) {
  ifelse(format(x, "%m") == "01" | (format(x, "%m") == "08" & format(x, "%Y") == "2021"), 
         format(x, "%b-%Y"), 
         format(x, "%b"))
}

co <- QAall %>%
  filter(PARAMETRO == "CO") |> 
  group_by(date) %>%
  summarise(VALOR = mean(VALOR))

ggplot(co, aes(x = date, y = VALOR)) +
  geom_line(color = "#2C3E50", size = 1) +  # Adjust line color and thickness
  scale_x_date(date_breaks = "1 month", labels = custom_date_labels2,
               expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  labs(
    title = "Daily Mean CO Level (2021-2024)", 
    y = "Concentration (mg/m³)"
  ) +
  theme_minimal(base_size = 12) +  # Use a minimal theme for a cleaner look
  theme(
    panel.grid.major.x = element_blank(),  # Remove vertical gridlines
    panel.grid.major.y = element_line(linewidth = 0.5, linetype = 'solid', colour = "lightgray"),  
    panel.grid.minor = element_blank(),  # Remove minor gridlines
    axis.line = element_line(colour = "black"),  # Keep axis lines
    axis.text.x = element_text(size = 11, angle = 45, hjust = 1),  # Rotate x-axis labels slightly
    axis.text.y = element_text(size = 11),
    axis.title = element_text(size = 13),
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),  # Emphasize title
    plot.subtitle = element_text(size = 14, hjust = 0.5),  # Add a subtitle
    plot.caption = element_text(size = 10, hjust = 1),  # Add a data source caption
    axis.title.x = element_blank()
  )
Figure 8. Daily mean CO level.

Figure 8. Daily mean CO level.