Extended vignette for EEAaq: Handle Air Quality Data from the European Environment Agency Data Portal
Extended vignette
Paolo Maranzano, University of Milano-Bicocca, Italy, paolo.maranzano@unimib.it
Riccardo Borgoni, University of Milano-Bicocca, Italy, riccardo.borgoni@unimib.it
Agostino Tassan Mazzocco, University of Milano-Bicocca, Italy
Samir Doghmi, University of Milano-Bicocca, Italy
20 gennaio 2026
The EEAaq package allows users to retrieve air quality data for multiple geographical zones, pollutants, and time periods in a single request. Queries are submitted as lists, which enables flexibility in specifying combinations of parameters.
library(EEAaq)
`%>%` <- dplyr::`%>%`EEAaq_get_data
Below we demonstrate the use of query by using different combinations of user-defined arguments.
Retrieve NO\(_2\) and PM\(_10\) data for the specified station IDs
#' ### Download PM10 data for the province (NUTS-3) of Milano (Italy) from January 1st 2024 to January 31st, 2025
IDstations <- EEAaq_get_stations(byStation = TRUE, complete = FALSE)
IDstations <- IDstations %>%
dplyr::filter(NUTS3 %in% c("Milano")) %>%
dplyr::pull(AirQualityStationEoICode) %>%
unique()
data <- EEAaq_get_data(IDstations = IDstations, pollutants = c("PM10", "NO2"),
from = "2024-01-01", to = "2025-01-31", verbose = TRUE)## Inputs check started at 2026-01-20 18:01:31.468185## Warning in FUN(X[[i]], ...): The following stations: IT0480A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT0592A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT0593A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT0623A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT0694A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT0761A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT1174A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT1203A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT1290A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT2232A are not associated
## with the pollutant : PM10## Warning in FUN(X[[i]], ...): The following stations: IT0524A are not associated
## with the pollutant : NO2## Download started at 2026-01-20 18:01:31.692346
## [[1]]
## [1] "{\"countries\": [\"IT\"],\"cities\": [],\"pollutants\": [\"PM10\", \"NO2\"],\"dataset\": \"2\",\"dateTimeStart\": \"2024-01-01T00:00:00Z\",\"dateTimeEnd\": \"2024-12-31T00:00:00Z\"}"
##
## [[2]]
## [1] "{\"countries\": [\"IT\"],\"cities\": [],\"pollutants\": [\"PM10\", \"NO2\"],\"dataset\": \"1\",\"dateTimeStart\": \"2025-01-01T00:00:00Z\",\"dateTimeEnd\": \"2025-01-31T00:00:00Z\"}"
##
## Post-processing started at 2026-01-20 18:02:05.734441## Warning in EEAaq_get_data(IDstations = IDstations, pollutants = c("PM10", : The
## following stations are not included in stations dataframe : IT0524A, IT0593A,
## IT0623A# Preview the first few rows of the dataset
head(data)unique(data$AirQualityStationEoICode)## [1] "IT0477A" "IT0480A" "IT0592A" "IT0694A" "IT0705A" "IT0706A" "IT0761A"
## [8] "IT1010A" "IT1016A" "IT1174A" "IT1203A" "IT1290A" "IT1294A" "IT1692A"
## [15] "IT1792A" "IT2232A"
Note 1: If the query’s IDstations parameter corresponds
to a valid CITY_NAME (i.e., not NULL in the dataset), the
function will return the corresponding data. If no valid CITY_NAME is
associated with the IDstations, the function attempts to retrieve all
available data for the entire country and subsequently filter for the
specified IDstations.
Note 3: If the parameters used in the query include
polygon or quadrant, the function outputs an
EEAaq_df_sfc object. Otherwise, it returns an EEAaq_df object, which is
a tibble dataframe.
EEAaq map stations
EEAaq_map_stations generates a static or dynamic map of
user-defined monitoring stations. The function accepts as input either
an object of the EEAaq_df class (default output of the
EEAaq_get_data function), or all other parameters
specifying the area and the pollutants.
Map the stations using as EEAaq_df object, which
contains NO\(_2\) PM\(_{10}\) data for Italy.
# Static map of available stations across the whole country. External borders are given by the
#' ### union of the available regions (NUTS-2), while municipalities (LAUs) are used as inner borders.
EEAaq_map_stations(data = data,
NUTS_extborder = "NUTS2", NUTS_intborder = "LAU",
color = TRUE, dynamic = FALSE)
Note: The external borders are given by the union of
the available regions (NUTS-2), while municipalities (LAUs) are used as
inner borders.
Dynamic (interactive) map of all the NO\(_2\) and PM\(_{10}\) monitoring stations in Milano
EEAaq_map_stations(data = data,
NUTS_extborder = "NUTS2", NUTS_intborder = "NUTS3",
color = TRUE, dynamic = TRUE)EEAaq summary
This function aims to describe the dataset that has been previously
imported, both at a global level, which means considering the complete
set of time stamps and monitoring stations in the dataset, and at the
station-specific level, where summary statistics and information are
grouped by monitoring station.
In addition to basic exploratory descriptive statistics (e.g., average
pollutant concentration, variability, measures of skewness and
kurtosis), the function provides information about the gap length and
the correlation between pollutants if at least two pollutants are
considered simultaneously.
The EEAaq_summary function receives as input an
EEAaq_df object, i.e. the output of the EEAaq get data
function.
Compute the descriptive statistics
summ <- EEAaq_summary(data = data)## The dataset contains:
## ** 155216 total observations
## ** 16 stations
## ** 9503 time stamps: from 2024-01-01 01:00:00 to 2025-01-31Print screen the global statsitics
summ$SummaryPrint screen the station-specific statsitics
summ$Summary_byStat$Mean_byStatPrint screen the linear correlation matrix
summ$Corr_MatrixEEAaq time aggregate
Recall that most pollutants are monitored by EEA on a hourly or daily
basis, posing challenges for interpretation and representation. The
EEAaq_time_aggregate function simplifies this by
aggregating data into annual, monthly, weekly, daily, or hourly
intervals, generating summary statistics for each station in an
EEAaq_taggr_df object.
Get the station-specific monthly minimum, maximum, average and median concentrations of NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands
t_aggr <- EEAaq_time_aggregate(
data = data,
frequency = "monthly",
aggr_fun = c("min", "max", "mean", "median" )
)Print screen of the aggregated (monthly) data
t_aggr$TimeAggrPrint screen of the PM\(_{10}\) aggregated data only
t_aggr$TimeAggr_byPollutant$PM10EEAaq_idw_map
To enable quick and intuitive visual analysis, the
EEAaq_idw_map function provides spatial interpolation maps
using the Inverse Distance Weighting (IDW) method (Shepard, 1968). This
technique estimates the value of a variable at unknown locations by
calculating a weighted average of known values, with weights inversely
proportional to the distance from known points. Closer points contribute
more heavily to the estimate, making it a practical approach for
interpolating geolocated air quality data.
Generate IDW interpolated maps of monthly average concentrations of NO\(_2\) in iTALY
EEAaq::EEAaq_idw_map(data = t_aggr, pollutant = "NO2", aggr_fun = "mean",
distinct = TRUE, gradient = FALSE,
dynamic = FALSE,
NUTS_filler = "NUTS3", NUTS_extborder = "NUTS2")## Map initialization started at 2026-01-20 18:02:46.626032
## Map initialization ended at 2026-01-20 18:02:56.953296
## Computing IDW interpolation started at 2026-01-20 18:02:56.953478
## Computing IDW interpolation for: 2024-01-01, 1 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-02-01, 2 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-03-01, 3 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-04-01, 4 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-05-01, 5 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-06-01, 6 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-07-01, 7 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-08-01, 8 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-09-01, 9 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-10-01, 10 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-11-01, 11 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-12-01, 12 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2025-01-01, 13 of 13
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2026-01-20 18:03:52.406892## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
t_aggr$TimeAggr_byPollutant$NO2Generate IDW interpolated maps of the maximum monthly concentrations of NO\(_2\) in january and february 2024 in Italy
EEAaq::EEAaq_idw_map(
data = t_aggr$TimeAggr_byPollutant$NO2 %>% dplyr::filter(Date %in% c("2024-01-01","2024-02-01")),
pollutant = "NO2",
aggr_fun = "max",
distinct = TRUE,
gradient = TRUE,
idp = 2,NUTS_extborder = "NUTS2",NUTS_filler = "LAU"
)## Map initialization started at 2026-01-20 18:04:00.370651
## Map initialization ended at 2026-01-20 18:04:11.973543
## Computing IDW interpolation started at 2026-01-20 18:04:11.973739
## Computing IDW interpolation for: 2024-01-01, 1 of 2
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2024-02-01, 2 of 2
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2026-01-20 18:04:24.030354## [[1]]
##
## [[2]]