10.2 Data!

10.2.1 Getting it

We get them from European Environment Agency

We are using annual averages of \(PM_2.5\) recorded at monitoring stations (2015, 2016,2017).

my_colnames <- c("ID", "CountryOrTerritory",    "ReportingYear",    "UpdateTime",
                 "StationLocalId",  "SamplingPointLocalId", "SamplingPoint_Latitude",
                 "SamplingPoint_Longitude", "Pollutant",    "AggregationType",
                 "Namespace",   "Unit", "BeginPosition",    "EndPosition",
                 "Validity",    "Verification", "DataCoverage", "DataCapture",
                 "TimeCoverage",    "AQValue")


d <- read.csv("https://www.paulamoraga.com/assets/datasets/dataPM25.csv", 
              skip = 1, 
              header = FALSE,
              col.names = my_colnames)

10.2.2 Cleaning it

Data are in 4326 (lat/long)

d <- d[, c(
  "ReportingYear", "StationLocalId",
  "SamplingPoint_Longitude",
  "SamplingPoint_Latitude",
  "AQValue"
)]
names(d) <- c("year", "id", "long", "lat", "value")

p <- st_as_sf(data.frame(long = d$long, lat = d$lat),
              coords = c("long", "lat"))
st_crs(p) <- st_crs(4326)
p <- st_transform(p, 25830)
d[, c("x", "y")] <- st_coordinates(p) # new coords in 25830

ind <- st_intersects(m, p) # return a sparse geometry binary predicate 
d <- d[ind[[1]], ]  # it is a list with just "spain" so we simplify [[]] it

We plot it!

ggplot(m) + geom_sf() + coord_sf(datum = NA) +
  geom_point(
    data = d, aes(x = x, y = y, color = value),
    size = 2
  ) +
  labs(x = "", y = "") +
  scale_color_viridis() +
  facet_wrap(~year) +
  theme_bw()