Get stations catalogue

The list of monitoring stations can be downloaded using the function ukair_catalogue() with no input parameters, as in the example below.

# Get full catalogue
stations <- ukair_catalogue()
stations
#> # A tibble: 6,608 x 16
#>    UK.AIR.ID EU.Site.ID EMEP.Site.ID Site.Name Environment.Type Zone 
#>    <chr>     <chr>      <chr>        <chr>     <chr>            <chr>
#>  1 UKA15910  <NA>       <NA>         SH15, SH… <NA>             Sout…
#>  2 UKA15956  <NA>       <NA>         SH15, SH… <NA>             Sout…
#>  3 UKA16663  <NA>       <NA>         SH15, SH… <NA>             Sout…
#>  4 UKA16097  <NA>       <NA>         10 High … <NA>             Sout…
#>  5 UKA12536  <NA>       <NA>         ABBOTS L… <NA>             Grea…
#>  6 UKA12949  <NA>       <NA>         ABERBARG… <NA>             Sout…
#>  7 UKA12399  <NA>       <NA>         ABERCARN… <NA>             Sout…
#>  8 UKA13340  <NA>       <NA>         ABERDARE… <NA>             Sout…
#>  9 UKA13341  <NA>       <NA>         ABERDARE… <NA>             Sout…
#> 10 UKA15369  <NA>       <NA>         ABERDARE… <NA>             Sout…
#> # … with 6,598 more rows, and 10 more variables: Start.Date <dttm>,
#> #   End.Date <dttm>, Latitude <dbl>, Longitude <dbl>, Northing <dbl>,
#> #   Easting <dbl>, Altitude..m. <dbl>, Networks <chr>,
#> #   AURN.Pollutants.Measured <chr>, Site.Description <chr>

There are currently 6608 stations in UK-AIR. The same function, can be used to filter the catalogue using the following input parameters:

• site_name IDs of specific site (UK.AIR.ID). By default this is left blank to get info on all the available sites.
• pollutant This is an integer between 1 and 10. Default is 9999, which means all the pollutants.
• group_id This is the identification number of a group of stations. Default is 9999 which means all available networks.
• closed This is set to TRUE to include closed stations, FALSE otherwise.
• country_id This is the identification number of the country, it can be an integer between 1 and 6. Default is 9999, which means all the countries.
• region_id This is the identification number of the region. 1 = Aberdeen City, etc. (for the full list see https://uk-air.defra.gov.uk/). Default is 9999, which means all the local authorities.

stations_EnglandOzone <- ukair_catalogue(pollutant = 1, country_id = 1)
stations_EnglandOzone
#> # A tibble: 108 x 16
#>    UK.AIR.ID EU.Site.ID EMEP.Site.ID Site.Name Environment.Type Zone 
#>    <chr>     <chr>      <chr>        <chr>     <chr>            <chr>
#>  1 UKA00353  GB0681A    <NA>         Barnsley… Background Urban York…
#>  2 UKA00626  GB1067A    <NA>         Birmingh… Traffic Urban    West…
#>  3 UKA00559  GB1013A    <NA>         Birmingh… Background Urban West…
#>  4 UKA00214  GB0569A    <NA>         Birmingh… Background Urban West…
#>  5 UKA00229  GB0595A    <NA>         Birmingh… Background Urban West…
#>  6 UKA00655  GB1097A    <NA>         Birmingh… Background Urban West…
#>  7 UKA00479  GB0851A    <NA>         Birmingh… Background Urban West…
#>  8 UKA00543  GB0960A    <NA>         Birmingh… Traffic Urban    West…
#>  9 UKA00412  GB0727A    <NA>         Blackpool Background Urban Blac…
#> 10 UKA00488  GB0882A    <NA>         Blackpoo… Background Urban Blac…
#> # … with 98 more rows, and 10 more variables: Start.Date <dttm>,
#> #   End.Date <dttm>, Latitude <dbl>, Longitude <dbl>, Northing <dbl>,
#> #   Easting <dbl>, Altitude..m. <dbl>, Networks <chr>,
#> #   AURN.Pollutants.Measured <chr>, Site.Description <chr>

The example above shows how to retrieve the 108 stations in England in which ozone is measured.

Get missing coordinates

Locating a station is extremely important to be able to carry out any spatial analysis. If coordinates are missing, for some stations in the catalogue, it might be possible to retrieve Easting and Northing (coordinates in the British National Grid) from DEFRA’s web pages, transform them to latitude and longitude and populate the missing coordinates as shown below.

# How many stations have missing coordinates?
length(which(is.na(stations$Latitude) | is.na(stations$Longitude)))
#> [1] 2763

# Scrape DEFRA website to get Easting/Northing (if available)
stations <- ukair_get_coordinates(stations)

# How many stations still have missing coordinates?
length(which(is.na(stations$Latitude) | is.na(stations$Longitude)))
#> [1] 2

Check hourly data availability

Pollution data started to be collected in 1972 and consists of hourly concentration of various species (in micrograms/m³), such as ozone (O₃), particulate matters (PM_2.5 and PM₁₀), nitrogen dioxide (NO₂), sulphur dioxide (SO₂), and so on.

The ID under which these data are available differs from the UK.AIR.ID. The catalogue does not contain this additional station ID (called SiteID hereafter) but DEFRA’s web pages contain references to both the UK.AIR.ID and the SiteID. The function below uses as input the UK.AIR.ID and outputs the SiteID, if available.

stations$SiteID <- ukair_get_site_id(stations$UK.AIR.ID)

Please note this function takes several minutes to run.

Cached catalogue

For convenience, a cached version of the catalogue (last updated in December 2019) is included in the package and can be loaded using the following command:

data("stations")

The cached catalogue contains all the available siteIDs and coordinates and can be used offline as lookup table to find out the correspondence between the UK.AIR.ID and SiteID, as well as to investigate station characteristics.

Get hourly data

Once the SiteID is known, time series for a given station can be retrieved in one line of code:

# Get 1 year of hourly ozone data from London Marylebone Road monitoring station
df <- ukair_get_hourly_data("MY1", years = 2015)

# Aggregate to daily means and plot
# please note we use the zoo package here because time series could be irregular
library("zoo")
my1 <- zoo(x = df$Ozone, order.by = as.POSIXlt(df$datetime))

daily_means <- aggregate(my1, as.Date(as.POSIXlt(df$datetime)), mean)

plot(daily_means, main = "", xlab = "",
     ylab = expression(paste("Ozone concentration [", mu, "g/", m^3, "]")))

Hourly ozone data from London Marylebone Road monitoring station in 2015

The above figure shows the highest concentrations happen in late spring and at the beginning of summer. In order to check whether this happens every year, we can download multiple years of data and then compare them.

The code below explores the distribution of ozone by month. The resulting box plots show that the highest concentrations usually occurr during April/May and that these vary year-by-year.

# Get 15 years of hourly ozone data from the same monitoring station
library("ggplot2")
library("dplyr")
library("lubridate")

df <- ukair_get_hourly_data("MY1", years = 2000:2015)
df <- mutate(df, 
             year = year(datetime),
             month = month(datetime),
             year_month = strftime(datetime, "%Y-%m"))

df %>%
  group_by(month, year_month) %>%
  summarize(ozone = mean(Ozone, na.rm=TRUE)) %>%
  ggplot() +
  geom_boxplot(aes(x = as.factor(month), y = ozone, group = month),
               outlier.shape = NA) +
  xlab("Month of the year") +
  ylab(expression(paste("Ozone concentration (", mu, "g/",m^3,")"))) +
  ggtitle("15 years of hourly ozone data from London Marylebone Road monitoring station")

Plotting stations’ locations

After scraping DEFRA’s web pages, almost all the stations have valid coordinates. You can create an interactive map using leaflet. The code below generates a map where blue circles are all the stations with valid coordinates, while red circles show locations with available hourly data.

# Keep only station with coordinates
stations_with_coords <- stations[complete.cases(stations[, c("Longitude",
                                                             "Latitude")]), ]
# Keep only station with known SiteID
stations_with_SiteID <- which(!is.na(stations_with_coords$SiteID))

# An interactive map
library("leaflet")
leaflet(data = stations_with_coords) %>% addTiles() %>% 
  addCircleMarkers(lng = ~Longitude, 
                   lat = ~Latitude,  
                   popup = ~SiteID,
                   radius = 1, color="blue", fill = FALSE) %>%
  addCircleMarkers(lng = ~Longitude[stations_with_SiteID], 
                   lat = ~Latitude[stations_with_SiteID], 
                   radius = 0.5, color="red", 
                   popup = ~SiteID[stations_with_SiteID])

Analyse the spatial distribution of the monitoring stations

Below are two plots showing the spatial distribution of the monitoring stations. These are concentrated largely in urban areas and mostly estimate the background level of concentration of pollutants.

# Zone
dotchart(as.matrix(table(stations$Zone))[,1])

Spatial distribution of the monitoring stations across zones.

# Environment.Type
dotchart(as.matrix(table(stations$Environment.Type[stations$Environment.Type != "Unknown Unknown"]))[,1])

Spatial distribution of the monitoring stations across environment types.

Use multiple cores to speed up data retrieval from numerous sites

The acquisition of data from hundreds of sites takes only few minutes:

library("parallel")
 
# Use detectCores() to find out many cores are available on your machine
cl <- makeCluster(getOption("cl.cores", detectCores()))

system.time(myList <- parLapply(cl, stations$SiteID[stations_with_SiteID], 
                                ukair_get_hourly_data, years=1999:2016))

stopCluster(cl)

df <- bind_rows(myList)