###### Paper: ARPALData: an R package for retrieving and analyzing air quality and weather data from ARPA Lombardia (Italy). ###### Authors: After the revision process ###### Journal: Environmental and Ecological Statistics ###### Date: 25/09/2023 ##### The following code reproduces the three guided examples discussed at Section 5 of the paper ##### Setup library(ARPALData) library(tidyverse) library(ggplot2) library(ggpubr) library(latex2exp) setwd("~/Ricerca/ARPALData_package") '%notin%' <- Negate('%in%') #################################################################### ########## Example 1: Characterizing AQ by local contexts ########## #################################################################### ### Download registry table of air quality monitoring sites reg <- get_ARPA_Lombardia_AQ_registry() ### Filter sensors: # Stations activated before 2019 # Currently on service / active # Measuring PM10 reg <- reg %>% filter(DateStart <= "2019-01-01", is.na(DateStop), Pollutant %in% c("PM10")) ### Stations type included in the sample table(reg$ARPA_stat_type) # Suburban-Traffic: 1 station # Suburban-Industrial: 3 stations # Urban-Industrial: 3 stations # Rural-Industrial: 2 stations ### Remove Suburban-Traffic, Suburban-Industrial, Urban-Industrial and Rural-Industrial due to low sample sizes reg <- reg %>% filter(ARPA_stat_type %notin% c("ST","SI","UI","RI")) ### Download weekly average and weekly maximum of PM10 for all the available stations from 2019 to 2022 # Parallel computation activated data <- get_ARPA_Lombardia_AQ_data(ID_station = reg$IDStation, Date_begin = "2019-01-01", Date_end = "2022-12-31", Frequency = "weekly", Var_vec = c("PM10","PM10"), Fns_vec = c("mean","max"), parallel = TRUE) ### Associating metadata (ARPA type classification and ARPA zoning) to observations reg_red <- reg %>% select(IDStation,ARPA_zone,ARPA_stat_type) data <- left_join(x = data, y = reg_red, by = c("IDStation")) data <- data %>% mutate(Date = lubridate::ymd(Date)) # save(data, reg_red, file="Data_ex1.RData") load("~/Ricerca/ARPALData_package/Data_ex1.RData") ### Compute weekly descriptive statistics for average PM10 concentrations by station type # Weekly minimum of average PM10 by station type # Weekly mean of average PM10 by station type # Weekly maximum of average PM10 by station type data_by_type <- data %>% group_by(Date,ARPA_stat_type) %>% summarise( min_avgPM10 = min(PM10_mean, na.rm = T), mean_avgPM10 = mean(PM10_mean, na.rm = T), max_avgPM10 = max(PM10_mean, na.rm = T) ) %>% ungroup() p1 <- data_by_type %>% pivot_longer(cols = 3:last_col(), names_to = "Pollutant", values_to = "Concentrations") %>% ggplot(mapping = aes(x = Date, col = Pollutant)) + geom_line(mapping = aes(y = Concentrations)) + facet_wrap(~ ARPA_stat_type) + labs(x = "", y = latex2exp::TeX("$\\mu g/m^3$"), title = latex2exp::TeX("Weekly average PM$_{10}$ concentrations ($\\mu g/m^3$) by stations type")) + theme_bw() ggpubr::ggexport(p1,width = 1800, height = 1200, res = 150, filename = "ARPALData_Ex1_ByType.png") ### Compute weekly descriptive statistics for maximum PM10 concentrations by zoning # Weekly minimum of maximum PM10 by zoning # Weekly mean of maximum PM10 by zoning # Weekly maximum of maximum PM10 by zoning zoning <- ARPALData::get_ARPA_Lombardia_zoning(plot_map = FALSE) data_by_zone <- data %>% group_by(Date,ARPA_zone) %>% summarise( min_maxPM10 = min(PM10_max, na.rm = T), mean_maxPM10 = mean(PM10_max, na.rm = T), max_maxPM10 = max(PM10_max, na.rm = T) ) %>% ungroup() data_by_zone <- left_join(x = data_by_zone, y = zoning %>% select(Cod_Zone,Zone), by = c("ARPA_zone" = "Cod_Zone")) p1 <- data_by_zone %>% pivot_longer(cols = contains("PM10"), names_to = "Pollutant", values_to = "Concentrations") %>% ggplot(mapping = aes(x = Date, col = Pollutant)) + geom_line(mapping = aes(y = Concentrations)) + facet_wrap(~ Zone) + labs(x = "", y = latex2exp::TeX("$\\mu g/m^3$"), title = latex2exp::TeX("Weekly maximum PM$_{10}$ concentrations ($\\mu g/m^3$) by ARPA zoning")) + theme_bw() + theme(legend.position = "") p2 <- ggplot(data = zoning, aes(fill = Zone)) + geom_sf() + labs(title = "ARPA Lombardia zoning", x = "Longitude", y = "Latitude") + theme(legend.position = "bottom", legend.text = element_text(size = 9)) + guides(fill=guide_legend(nrow=7,byrow=TRUE)) # Combining plots p_comb <- ggpubr::ggarrange(p1,p2, ncol = 2, nrow = 1, widths = c(2,1)) print(p_comb) ggpubr::ggexport(p_comb,width = 1800, height = 1200, res = 150, filename = "ARPALData_Ex1_ByZone.png") ############################################################################### ########## Example 2: AQ during COVID-19 lockdown at municipal level ########## ############################################################################### ### Download registry table of air quality monitoring sites reg_AQ <- get_ARPA_Lombardia_AQ_municipal_registry() ### Filter sensors: # Stations active in the region (currently on service) # Measuring NO2 concentrations reg_AQ <- reg_AQ %>% filter(is.na(DateStop), Pollutant %in% c("NO2")) ### Download daily NO2 concentrations at municipal level from 2018 to 2021 Data_mun <- get_ARPA_Lombardia_AQ_municipal_data( Date_begin = "2018-01-01",Date_end = "2021-12-31", Frequency = "daily", Var_vec = c("NO2_mean"), Fns_vec = c("mean"), verbose = TRUE ) # save(Data_mun, reg_AQ, file="Data_ex2.RData") load("~/Ricerca/ARPALData_package/Data_ex2.RData") ### Computing period averages (8th March - 18th May) of NO2 concentrations from 2018 to 2021 Data_spring <- Data_mun %>% filter(Date >= "2021-03-08" & Date <= "2021-05-18" | Date >= "2020-03-08" & Date <= "2020-05-18" | Date >= "2019-03-08" & Date <= "2019-05-18" | Date >= "2018-03-08" & Date <= "2018-05-18") Data_y <- Time_aggregate( Dataset = Data_spring, Frequency = "yearly" ) ### Compute reference value for the mean: average NO2 concentrations in 2018 mid_conc_2018 <- Data_y %>% filter(lubridate::year(Date) == 2018) %>% summarise(mean(NO2_mean,na.rm = T)) %>% pull() mid_conc_2018 ### Generate maps of average NO2 concentrations during the subperiod (3rd March - 8th May) from 2018 to 2021 map_18 <- ARPALdf_Summary_map( Data = Data_y %>% filter(lubridate::year(Date) == 2018), Title_main = latex2exp::TeX("March 8$^{th}$ to May 18$^{th}$, 2018"), Variable = "NO2_mean", val_midpoint = mid_conc_2018 ) map_19 <- ARPALdf_Summary_map( Data = Data_y %>% filter(lubridate::year(Date) == 2019), Title_main = latex2exp::TeX("March 8$^{th}$ to May 18$^{th}$, 2019"), Variable = "NO2_mean", val_midpoint = mid_conc_2018 ) map_20 <- ARPALdf_Summary_map( Data = Data_y %>% filter(lubridate::year(Date) == 2020), Title_main = latex2exp::TeX("March 8$^{th}$ to May 18$^{th}$, 2020"), Variable = "NO2_mean", val_midpoint = mid_conc_2018 ) ### Map for 2021 map_21 <- ARPALdf_Summary_map( Data = Data_y %>% filter(lubridate::year(Date) == 2021), Title_main = latex2exp::TeX("March 8$^{th}$ to May 18$^{th}$, 2021"), Variable = "NO2_mean", val_midpoint = mid_conc_2018 ) fig_comb <- ggarrange(map_18,NULL,map_19,map_20,NULL,map_21, ncol = 3, nrow = 2, widths = c(1, 0.0005, 1, 1, 0.0005, 1), common.legend = T, legend = "bottom") fig_comb <- annotate_figure(p = fig_comb, top = text_grob(latex2exp::TeX("Average NO$_{2}$ concentrations by municipality during spring time"), col="blue",face = "bold",size = 14)) ggpubr::ggexport(fig_comb,width = 1800, height = 1200, res = 150, filename = "ARPALData_Ex2_Maps.png") ############################################################################################# ########## Example 3: characterization of meteorological phenomena across Lombardy ########## ############################################################################################# ### Downloading weather measurements for 2021 at monthly frequency # Cumulated rainfall, average temperature and maximum temperature data <- get_ARPA_Lombardia_W_data( ID_station = NULL, Date_begin = "2021-01-01",Date_end = "2021-12-31", Frequency = "monthly", Var_vec = c("Rainfall","Temperature","Wind_speed","Wind_direction"), Fns_vec = c("sum","mean","mean","mean"), parallel = TRUE ) # save(data, file="Data_ex3.RData") load("~/Ricerca/ARPALData_package/Data_ex3.RData") ### Summary statistics and plots for the dataset: # Descriptive statistics for the whole sample # Descriptive statistics for each station and for each year # Gap length (missing data) analysis # Outlier analysis # Correlation analysis for each station # Density plots for each variable summ_data <- ARPALdf_Summary( Data = data, by_IDStat = TRUE, by_Year = TRUE, gap_length = TRUE, outlier = TRUE, correlation = TRUE, density = TRUE, histogram = FALSE, verbose = FALSE ) library(xtable) summ_data$Hampel %>% View() head(summ_data$Gap_length$Rainfall,5) %>% View() reg_tex <- xtable(summ_data$Hampel, caption = "Outlier analysis for each variable (Hampel filter)", align=c("l","l","c","c","c","c","c","c")) head(summ_data$Gap_length$Wind_speed,100) %>% View() sum_stats_tex <- xtable(head(summ_data$Gap_length$Wind_speed,10), caption = "Gap length (missing data) analysis for 5 stations", align=c("l","l","l","c","c","c","c","c","c","c","c","c","c")) ### Subsampling the observations according to the climatic season: # spring/summer data (april to september) and autumn/winter data (october to march) data_spring_summer <- data %>% filter(lubridate::month(Date) %in% c(4,5,6,7,8,9)) data_autumn_winter <- data %>% filter(lubridate::month(Date) %in% c(10,11,12,1,2,3)) ### Summary statistics afor each subsample summ_data_spring_summer <- ARPALdf_Summary( data_spring_summer, gap_length = F, outlier = F,correlation = T,density = F, verbose = F) summ_data_autumn_winter <- ARPALdf_Summary( data_autumn_winter, gap_length = F, outlier = F,correlation = T,density = F, verbose = F) ### Plotting average temperature by station and season # Computing average temperature in 2021 as central reference value mean_temp <- mean(data$Temperature,na.rm = T) # Map for spring/summer period map_T_summer <- ARPALdf_Summary_map(summ_data_spring_summer$Descr_by_IDStat$Mean_by_stat, Variable = "Temperature", Title_main = "Average temperature in Spring/Summmer 2021", col_scale = c("#FFFF00","#FF9933","#FF0000"), val_midpoint = mean_temp) # Map for autumn/winter period map_T_winter <- ARPALdf_Summary_map(summ_data_autumn_winter$Descr_by_IDStat$Mean_by_stat, Variable = "Temperature", Title_main = "Average temperature in Autumn/Winter 2021", col_scale = c("#FFFF00","#FF9933","#FF0000"), val_midpoint = mean_temp) # Combining plots fig_comb <- ggarrange(map_T_summer,map_T_winter,ncol = 2, nrow = 1) ggpubr::ggexport(fig_comb,width = 1800, height = 1200, res = 150, filename = "ARPALData_Ex3_Temp.png") ### Plotting linear correlation among maximum temperature and rainfall by station and season # Map for spring/summer period map_corr1 <- ARPALdf_Summary_map(summ_data_spring_summer$Cor_matrix, Variable = "Temperature_Rainfall", Title_main = "Correlation among max temperature and rainfall in Spring/Summmer 2020", val_midpoint = 0) # Map for autumn/winter period map_corr2 <- ARPALdf_Summary_map(summ_data_autumn_winter$Cor_matrix, Variable = "Temperature_Rainfall", Title_main = "Correlation among max temperature and rainfall in Autumn/Winter 2020", val_midpoint = 0) # Combining plots fig_comb <- ggarrange(map_corr1,map_corr2,ncol = 2, nrow = 1) ggpubr::ggexport(fig_comb,width = 1800, height = 1200, res = 150, filename = "ARPALData_Ex3_Corr.png")