Get the ward boundaries

ukgrid      <- "+init=epsg:27700"
latlong     <- "+init=epsg:4326"

authorities <- data.frame(authority_name = c('Waltham Forest'),
                          stringsAsFactors = F)

###  getting a geojson of UK wards from governmant data portal

url                 <- 'https://opendata.arcgis.com/datasets/d5c9c1d89a5a44e9a7f88f182ffe5ba2_2.geojson'
wards               <- readOGR(dsn = url, layer = "d5c9c1d89a5a44e9a7f88f182ffe5ba2_2")
## OGR data source with driver: GeoJSON 
## Source: "https://opendata.arcgis.com/datasets/d5c9c1d89a5a44e9a7f88f182ffe5ba2_2.geojson", layer: "d5c9c1d89a5a44e9a7f88f182ffe5ba2_2"
## with 9130 features
## It has 12 fields
wards               <- spTransform(wards, ukgrid)
wards               <- wards[wards$lad16nm %in% authorities$authority_name,]

Import the air quality files

no2_2013_raster <- raster('T:/Projects/WalthamForest/Modelling/PostLAEI2013_WF/2013/ASCII/PostLAEI_2013_WFSpeedUpdate_NO2.asc')
pm25_2013_raster <- raster('T:/Projects/WalthamForest/Modelling/PostLAEI2013_WF/2013/ASCII/PostLAEI_2013_WFSpeedUpdate_PM25.asc')

no2_2020_raster <- raster('T:/Projects/WalthamForest/Modelling/PostLAEI2013_WF/2020/ASCII/no2_2020.asc')
pm25_2020_raster <- raster('T:/Projects/WalthamForest/Modelling/PostLAEI2013_WF/2020/ASCII/pm25_2020.asc')

Check rasters

plot(no2_2013_raster)
plot(wards, add=T)
title('2013 NO2')

plot(pm25_2013_raster)
plot(wards, add=T)
title('2013 PM25')

plot(no2_2020_raster)
plot(wards, add=T)
title('2020 NO2')

plot(pm25_2020_raster)
plot(wards, add=T)
title('2020 PM25')

Crop the rasters and check they look ok

wards_extent      <- extent(st_bbox(wards)$xmin-40, st_bbox(wards)$xmax+40, 
                      st_bbox(wards)$ymin-40, st_bbox(wards)$ymax+40)

no2_2013_raster   <- crop(no2_2013_raster, wards_extent)
pm25_2013_raster  <- crop(pm25_2013_raster, wards_extent)
no2_2020_raster   <- crop(no2_2020_raster, wards_extent)
pm25_2020_raster  <- crop(pm25_2020_raster, wards_extent)

plot(no2_2013_raster)
plot(wards, add=T)
title('2013 NO2')

plot(pm25_2013_raster)
plot(wards, add=T)
title('2013 PM25')

plot(no2_2020_raster)
plot(wards, add=T)
title('2020 NO2')

plot(pm25_2020_raster)
plot(wards, add=T)
title('2020 PM25')

Extract the concentrations per ward

wards$mean_no2_2013    <- extract(no2_2013_raster, wards, fun=mean)
wards$mean_pm25_2013   <- extract(pm25_2013_raster, wards, fun=mean)

wards$mean_no2_2020    <- extract(no2_2020_raster, wards, fun=mean)
wards$mean_pm25_2020   <- extract(pm25_2020_raster, wards, fun=mean)

Prep work for pltting

wards             <- st_as_sf(wards)

wards$wd16nm       <- as.character(wards$wd16nm)
wards$id          <- 1:nrow(wards)
ward_centroids    <- data.frame(id = as.character(wards$id),
                                  name = as.character(wards$wd16nm),
                                  st_coordinates(st_centroid(wards)))
## Warning in st_centroid.sf(wards): st_centroid assumes attributes are
## constant over geometries of x
names(ward_centroids)[3:4] <- c('x','y')

The wards

ggplot(wards, aes(colour = wd16nm)) +
  geom_sf() +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        legend.background = element_rect(size=0.2, linetype="solid", colour ="black")) +
  scale_colour_manual(values = rep('black', 20), name = 'Ward name', labels = paste(1:20, wards$wd16nm)) +
  geom_text(data = ward_centroids, aes(x, y, label = id), colour='black')

Mean 2013 NO2 concentrations

ggplot(wards) +
  geom_sf(aes(fill=mean_no2_2013)) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        legend.background = element_rect(size=0.2, linetype="solid", colour ="black"))

Mean 2020 PM2.5 concentrations

ggplot(wards) +
  geom_sf(aes(fill=mean_pm25_2013)) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        legend.background = element_rect(size=0.2, linetype="solid", colour ="black"))

Mean 2013 NO2 concentrations

ggplot(wards) +
  geom_sf(aes(fill=mean_no2_2020)) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        legend.background = element_rect(size=0.2, linetype="solid", colour ="black"))

Mean 2020 PM2.5 concentrations

ggplot(wards) +
  geom_sf(aes(fill=mean_pm25_2020)) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        legend.background = element_rect(size=0.2, linetype="solid", colour ="black"))

Write concentrations to a CSV

st_geometry(wards) <- NULL
write.csv(data.frame(wards[,c('wd16cd', 'wd16nm', 'mean_no2_2013', 'mean_pm25_2013', 'mean_no2_2020', 'mean_pm25_2020')]), 'wf_ward_concs.csv', row.names = F)

Download final CSV from here