R script to prepare air quality long term averages for Italian provinces

get_aqdata.R

# packages
library(saqgetr)
library(dplyr)
library(rgdal)
library(tidyr)
library(rgeos)
library(readr)

# extraction
data_annual <- get_saq_simple_summaries(summary = "annual_means") 
data_sites <- get_saq_sites() 

# selection of data
data_annual %>%
  filter(variable %in% c("no2", "pm10"#, "pm2.5"
                         )) %>%
  mutate(ndays=as.numeric(date_end - date),
         stepinaday=if_else(summary_source==20,1,24),
         nstep=ndays*stepinaday) %>%
  filter(count>=0.9*nstep)-> dat

# selection of sites
data_sites %>%
  dplyr::filter(country_iso_code=="IT",
         !is.na(site_type),
         !is.na(site_area),
         !is.na(date_start),
         elevation>-10,
         elevation<1000,
         site_type!="industrial",
         site_type!="traffic",     # uncomment this line if you need averages on kerbside stations as well
         substr(site_area,1,5)!="rural") -> sites

# add Italian provinces
prov <- readOGR("dat/province.geojson")
pr <- gBuffer(gSimplify(prov,0.01,topologyPreserve = T),byid = T,width = 0.01)   
prov <- SpatialPolygonsDataFrame(Sr = pr, data = prov@data, match.ID = T)
pts <- data.frame(x=sites$longitude, y=sites$latitude)
coordinates(pts) <- ~ x+y
pts@proj4string <- prov@proj4string
pp <- over(pts, prov)
sites$prov <- pp$NOME_PRO
sites %>%
  filter(!is.na(prov)) -> sites

# merging data with sites
left_join(sites %>% 
            dplyr::select(site:elevation, site_type, site_area, prov),
          dat %>%
            dplyr::select(date, site, variable, value)) -> Dat

# selecting period
firstyear <- 2013
lastyear <- 2018
nyrs <- lastyear-firstyear+1
Dat %>%
  mutate(year=as.numeric(format(date,"%Y"))) %>%
  filter(year>=firstyear, year<=lastyear) %>%
  group_by(variable, site) %>%
  filter(n()>=2) -> Dat

# spatial aggregation 
Dat %>%
  group_by(site_type, prov, year, variable) %>%
  summarize(value=mean(value)) -> Dat

# temporal aggregation
Dat %>%
  group_by(site_type, prov, variable) %>%
  filter(n()>=nyrs) %>%
  summarize(value=round(mean(value),2)) -> avDat

# table (annual averaged)
Dat %>%
  ungroup() %>%
  mutate(variable=paste(variable,site_type,sep="_"),
         value=round(value,2),
         site_type=NULL) %>%
  spread(variable,value) -> Tab
left_join(data.frame(prov=unique(prov@data$NOME_PRO)),
          Tab) %>%
  arrange(prov) %>%
  filter(!is.na(year)) -> Tab
write_csv(Tab,paste0("AQdata_",firstyear,"-",lastyear,"_annualAve.csv"))

# table (multiannual averaged)
avDat %>%
  ungroup() %>%
  mutate(variable=paste(variable,site_type,sep="_"),
         site_type=NULL) %>%
  spread(variable,value) -> Tab
left_join(data.frame(prov=unique(prov@data$NOME_PRO)),
          Tab) %>%
  arrange(prov)-> Tab
write_csv(Tab,paste0("AQdata_",firstyear,"-",lastyear,"_ave.csv"))