Waltham Forest Concentration Maps

Libraries

rm(list=ls())

### Loading any missing libraries

library(rvest, quietly = T)

## Warning: package 'rvest' was built under R version 3.5.1

## Warning: package 'xml2' was built under R version 3.5.1

library(stringr, quietly = T)

## Warning: package 'stringr' was built under R version 3.5.2

library(raster, quietly = T)

## Warning: package 'raster' was built under R version 3.5.3

## Warning: package 'sp' was built under R version 3.5.1

library(rgdal, quietly = T)

## Warning: package 'rgdal' was built under R version 3.5.1

## rgdal: version: 1.3-3, (SVN revision 759)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Program Files/R/R-3.5.0/library/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Program Files/R/R-3.5.0/library/rgdal/proj
##  Linking to sp version: 1.3-1

library(rgeos, quietly = T)

## Warning: package 'rgeos' was built under R version 3.5.2

## rgeos version: 0.4-2, (SVN revision 581)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE

library(ggplot2, quietly = T)

## Warning: package 'ggplot2' was built under R version 3.5.2

library(sf, quietly = T)

## Warning: package 'sf' was built under R version 3.5.2

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(rmarkdown, quietly = T)
library(rasterVis, quietly = T)

## Warning: package 'rasterVis' was built under R version 3.5.2

## Warning: package 'latticeExtra' was built under R version 3.5.1

## 
## Attaching package: 'latticeExtra'

## The following object is masked from 'package:ggplot2':
## 
##     layer

Data

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

wards               <-st_read('https://opendata.arcgis.com/datasets/d5c9c1d89a5a44e9a7f88f182ffe5ba2_2.geojson', quiet = TRUE)
wards               <- st_transform(wards, 27700)
wards               <- wards[wards$lad16nm %in% authorities$authority_name,]

Import the air quality files

no2_2013_raster  <- raster('../air_quality_inputs/no2_2013.asc')
nox_2013_raster  <- raster('../air_quality_inputs/nox_2013.asc')
pm25_2013_raster <- raster('../air_quality_inputs/pm25_2013.asc')
pm10_2013_raster <- raster('../air_quality_inputs/pm10_2013.asc')
  
no2_2020_raster  <- raster('../air_quality_inputs/no2_2020.asc')
nox_2020_raster  <- raster('../air_quality_inputs/nox_2020.asc')
pm10_2020_raster <- raster('../air_quality_inputs/pm10_2020.asc')
pm25_2020_raster <- raster('../air_quality_inputs/pm25_2020.asc')

Now import the ones that need cropping and re-gridding

no2_2013_minus_no_sc_raster       <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/no22013mnosc/w001001.adf')
nox_2013_minus_no_sc_raster       <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/nox2013mnosc/w001001.adf')
pm25_2013_minus_no_sc_raster      <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/pm252013mnosc/w001001.adf')
pm10_2013_minus_no_sc_raster      <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/pm102013mnosc/w001001.adf')

crs(no2_2013_minus_no_sc_raster)  <- CRS(ukgrid)
crs(nox_2013_minus_no_sc_raster)  <- CRS(ukgrid)
crs(pm25_2013_minus_no_sc_raster) <- CRS(ukgrid)
crs(pm10_2013_minus_no_sc_raster) <- CRS(ukgrid)

no2_2020_sc3_minus_sc1_raster  <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/no220sc3msc1/w001001.adf')
nox_2020_sc3_minus_sc1_raster  <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/nox20sc3msc1/w001001.adf')
pm25_2020_sc3_minus_sc1_raster <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/pm2520sc3msc1/w001001.adf')
pm10_2020_sc3_minus_sc1_raster <- raster('Z:/Projects/AutomatedAirPollutionMapping/CHASE/R_Mapping/Output/WF/Raster/pm1020sc3msc1/w001001.adf')

crs(no2_2020_sc3_minus_sc1_raster)  <- CRS(ukgrid)
crs(nox_2020_sc3_minus_sc1_raster)  <- CRS(ukgrid)
crs(pm25_2020_sc3_minus_sc1_raster) <- CRS(ukgrid)
crs(pm10_2020_sc3_minus_sc1_raster) <- CRS(ukgrid)

no2_2013_minus_no_sc_raster    <- crop(no2_2013_minus_no_sc_raster,    wards)
nox_2013_minus_no_sc_raster    <- crop(nox_2013_minus_no_sc_raster,    wards)
pm25_2013_minus_no_sc_raster   <- crop(pm25_2013_minus_no_sc_raster,   wards)
pm10_2013_minus_no_sc_raster   <- crop(pm10_2013_minus_no_sc_raster,   wards)

no2_2020_sc3_minus_sc1_raster  <- crop(no2_2020_sc3_minus_sc1_raster,  wards)
nox_2020_sc3_minus_sc1_raster  <- crop(nox_2020_sc3_minus_sc1_raster,  wards)
pm25_2020_sc3_minus_sc1_raster <- crop(pm25_2020_sc3_minus_sc1_raster, wards)
pm10_2020_sc3_minus_sc1_raster <- crop(pm10_2020_sc3_minus_sc1_raster, wards)

no2_2013_minus_no_sc_raster    <- mask(no2_2013_minus_no_sc_raster,    wards)
nox_2013_minus_no_sc_raster    <- mask(nox_2013_minus_no_sc_raster,    wards)
pm25_2013_minus_no_sc_raster   <- mask(pm25_2013_minus_no_sc_raster,   wards)
pm10_2013_minus_no_sc_raster   <- mask(pm10_2013_minus_no_sc_raster,   wards)

no2_2020_sc3_minus_sc1_raster  <- mask(no2_2020_sc3_minus_sc1_raster,  wards)
nox_2020_sc3_minus_sc1_raster  <- mask(nox_2020_sc3_minus_sc1_raster,  wards)
pm25_2020_sc3_minus_sc1_raster <- mask(pm25_2020_sc3_minus_sc1_raster, wards)
pm10_2020_sc3_minus_sc1_raster <- mask(pm10_2020_sc3_minus_sc1_raster, wards)

no2_2013_minus_no_sc_raster    <- disaggregate(no2_2013_minus_no_sc_raster,    fact = 4, method="bilinear")
nox_2013_minus_no_sc_raster    <- disaggregate(nox_2013_minus_no_sc_raster,    fact = 4, method="bilinear")
pm25_2013_minus_no_sc_raster   <- disaggregate(pm25_2013_minus_no_sc_raster,   fact = 4, method="bilinear")
pm10_2013_minus_no_sc_raster   <- disaggregate(pm10_2013_minus_no_sc_raster,   fact = 4, method="bilinear")

no2_2020_sc3_minus_sc1_raster  <- disaggregate(no2_2020_sc3_minus_sc1_raster,  fact = 4, method="bilinear")
nox_2020_sc3_minus_sc1_raster  <- disaggregate(nox_2020_sc3_minus_sc1_raster,  fact = 4, method="bilinear")
pm25_2020_sc3_minus_sc1_raster <- disaggregate(pm25_2020_sc3_minus_sc1_raster, fact = 4, method="bilinear")
pm10_2020_sc3_minus_sc1_raster <- disaggregate(pm10_2020_sc3_minus_sc1_raster, fact = 4, method="bilinear")

rm(wards, authorities)

Now import some colour files

source('https://raw.githubusercontent.com/KCL-ERG/colour_schemes/master/2013/pm25_laei2013_colours_breaks.R')
source('https://raw.githubusercontent.com/KCL-ERG/colour_schemes/master/2013/no2_laei2013_colours_breaks.R')
source('https://raw.githubusercontent.com/KCL-ERG/colour_schemes/master/2013/pm10_laei2013_colours_breaks.R')

Maps

NO2 2013

png('../map_outputs/no2_2013_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(no2_2013_raster,
          maxpixels = no2_2013_raster@ncols/2 * no2_2013_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(no2_2013_raster,plot)

NO2 2013

png('../map_outputs/no2_2020_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(no2_2020_raster,
          maxpixels = no2_2020_raster@ncols/2 * no2_2020_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(no2_2020_raster, plot)

NOx 2013

png('../map_outputs/nox_2013_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(nox_2013_raster,
          maxpixels = nox_2020_raster@ncols/2 * nox_2020_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(nox_2013_raster, plot)

NOx 2020

png('../map_outputs/nox_2020_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(nox_2020_raster,
          maxpixels = nox_2020_raster@ncols/2 * nox_2020_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = 17), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(nox_2020_raster, plot)

PM25 2013

png('../map_outputs/pm25_2013_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm25_2013_raster,
          maxpixels = pm25_2013_raster@ncols * pm25_2013_raster@nrows,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = 12),
            space = 'right',
            labels = list(at=seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = 12), 
                          labels = paste(" \n \n ",pm25_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm25_laei2013_colours,
          at = pm25_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm25_2013_raster, plot)

PM25 2020

png('../map_outputs/pm25_2020_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm25_2020_raster,
          maxpixels = pm25_2020_raster@ncols * pm25_2020_raster@nrows,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = 12),
            space = 'right',
            labels = list(at=seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = 12), 
                          labels = paste(" \n \n ",pm25_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm25_laei2013_colours,
          at = pm25_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm25_2020_raster, plot)

PM10 2013

png('../map_outputs/pm10_2013_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm10_2013_raster,
          maxpixels = pm10_2013_raster@ncols * pm10_2013_raster@nrows,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = 17),
            space = 'right',
            labels = list(at=seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = 17), 
                          labels = paste(" \n \n ",pm10_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm10_laei2013_colours,
          at = pm10_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm10_2013_raster, plot)

PM10 2020

png('../map_outputs/pm10_2020_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm10_2020_raster,
          maxpixels = pm10_2020_raster@ncols * pm10_2020_raster@nrows,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = 17),
            space = 'right',
            labels = list(at=seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = 17), 
                          labels = paste(" \n \n ",pm10_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm10_laei2013_colours,
          at = pm10_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm10_2020_raster, plot)

no2_2013_minus_no_sc_raster    <- no2_2013_minus_no_sc_raster    * -1
nox_2013_minus_no_sc_raster    <- nox_2013_minus_no_sc_raster    * -1
pm25_2013_minus_no_sc_raster   <- pm25_2013_minus_no_sc_raster   * -1
pm10_2013_minus_no_sc_raster   <- pm10_2013_minus_no_sc_raster   * -1

NO2 2013 minus no school run

no2_laei2013_breaks  <- c(0,format(round(quantile(no2_2013_minus_no_sc_raster, seq(0,1,length.out = 15)),4), scientific=F), 1) #17

no2_laei2013_labels  <- c("", paste(no2_laei2013_breaks[1:length(no2_laei2013_breaks)-1], "-", 
                                    no2_laei2013_breaks[2:length(no2_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/no2_2013_minus_no_school_run_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(no2_2013_minus_no_sc_raster,
          maxpixels = no2_2013_minus_no_sc_raster@ncols/2 * no2_2013_minus_no_sc_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(no2_2013_minus_no_sc_raster, plot)

NO2 2020 scenario 3 minus scenario 1

no2_laei2013_breaks  <- c(-1,format(round(quantile(no2_2020_sc3_minus_sc1_raster, seq(0,1,length.out = 15)),3), scientific=F), 0) #17

#no2_laei2013_colours <- c("#064AF4", "#0C95E9", "#19CFD2", "#82FDCF", "#68DE85", "#A4EB50", "#FFFF80", "#FFD600", "#FFAD5B", "#F97C00") #16 (one less above)

no2_laei2013_labels  <- c("", paste(no2_laei2013_breaks[1:length(no2_laei2013_breaks)-1], "to", 
                                    no2_laei2013_breaks[2:length(no2_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/no2_2020_sc3_minus_sc1_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(no2_2020_sc3_minus_sc1_raster,
          maxpixels = no2_2020_sc3_minus_sc1_raster@ncols/2 * no2_2020_sc3_minus_sc1_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(no2_2020_sc3_minus_sc1_raster, plot)

NOX 2013 minus no school run

no2_laei2013_breaks  <- c(0,format(round(quantile(nox_2013_minus_no_sc_raster, seq(0,1,length.out = 15)),4), scientific=F), 4) #17

#no2_laei2013_colours <- c("#064AF4", "#0C95E9", "#19CFD2", "#82FDCF", "#68DE85", "#A4EB50", "#FFFF80", "#FFD600", #"#FFAD5B", "#F97C00") #16 (one less above)

no2_laei2013_labels  <- c("", paste(no2_laei2013_breaks[1:length(no2_laei2013_breaks)-1], "-", 
                                    no2_laei2013_breaks[2:length(no2_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/nox_2013_minus_no_school_run_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(nox_2013_minus_no_sc_raster,
          maxpixels = nox_2013_minus_no_sc_raster@ncols/2 * nox_2013_minus_no_sc_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(nox_2013_minus_no_sc_raster, plot)

NOX 2020 scenario 3 minus scenario 1

no2_laei2013_breaks  <- c(-2,format(round(quantile(nox_2020_sc3_minus_sc1_raster, seq(0,1,length.out = 15)),3), scientific=F), 0) #17

no2_laei2013_labels  <- c("", paste(no2_laei2013_breaks[1:length(no2_laei2013_breaks)-1], "to", 
                                    no2_laei2013_breaks[2:length(no2_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/nox_2020_sc3_minus_sc1_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(nox_2020_sc3_minus_sc1_raster,
          maxpixels = nox_2020_sc3_minus_sc1_raster@ncols/2 * nox_2020_sc3_minus_sc1_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(no2_laei2013_breaks), max(no2_laei2013_breaks), length = length(no2_laei2013_breaks)), 
                          labels = paste(" \n \n ",no2_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = no2_laei2013_colours,
          at = no2_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(nox_2020_sc3_minus_sc1_raster, plot)

PM25 2013 minus no school run

pm25_laei2013_breaks  <- c(0,format(round(quantile(pm25_2013_minus_no_sc_raster, seq(0,1,length.out = 10)),4), scientific=F), 0.13) #17

pm25_laei2013_labels  <- c("", paste(pm25_laei2013_breaks[1:length(pm25_laei2013_breaks)-1], "-", 
                                    pm25_laei2013_breaks[2:length(pm25_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/pm25_2013_minus_no_school_run_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm25_2013_minus_no_sc_raster,
          maxpixels = pm25_2013_minus_no_sc_raster@ncols/2 * pm25_2013_minus_no_sc_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = length(pm25_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = length(pm25_laei2013_breaks)), 
                          labels = paste(" \n \n ",pm25_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm25_laei2013_colours,
          at = pm25_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm25_2013_minus_no_sc_raster, plot)

PM25 2020 scenario 3 minus scenario 1

pm25_laei2013_breaks  <- c(-0.06,format(round(quantile(pm25_2020_sc3_minus_sc1_raster, seq(0,1,length.out = 10)),4), scientific=F), 0) #17

pm25_laei2013_labels  <- c("", paste(pm25_laei2013_breaks[1:length(pm25_laei2013_breaks)-1], "to", 
                                    pm25_laei2013_breaks[2:length(pm25_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/pm25_2020_sc3_minus_sc1_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm25_2020_sc3_minus_sc1_raster,
          maxpixels = pm25_2020_sc3_minus_sc1_raster@ncols/2 * pm25_2020_sc3_minus_sc1_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = length(pm25_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(pm25_laei2013_breaks), max(pm25_laei2013_breaks), length = length(pm25_laei2013_breaks)), 
                          labels = paste(" \n \n ",pm25_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm25_laei2013_colours,
          at = pm25_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm25_2020_sc3_minus_sc1_raster, plot)

PM10 2013 minus no school run

pm10_laei2013_breaks  <- c(0,format(round(quantile(pm10_2013_minus_no_sc_raster, seq(0,1,length.out = 15)),4), scientific=F), 0.3) #17

pm10_laei2013_labels  <- c("", paste(pm10_laei2013_breaks[1:length(pm10_laei2013_breaks)-1], "-", 
                                    pm10_laei2013_breaks[2:length(pm10_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/pm10_2013_minus_no_school_run_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm10_2013_minus_no_sc_raster,
          maxpixels = pm10_2013_minus_no_sc_raster@ncols/2 * pm10_2013_minus_no_sc_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = length(pm10_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = length(pm10_laei2013_breaks)), 
                          labels = paste(" \n \n ",pm10_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm10_laei2013_colours,
          at = pm10_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm10_2013_minus_no_sc_raster, plot)

PM10 2020 scenario 3 minus scenario 1

pm10_laei2013_breaks  <- c(-0.014,format(round(quantile(pm10_2020_sc3_minus_sc1_raster, seq(0,1,length.out = 10)),4), scientific=F), 0) #17

pm10_laei2013_labels  <- c("", paste(pm25_laei2013_breaks[1:length(pm25_laei2013_breaks)-1], "to", 
                                    pm25_laei2013_breaks[2:length(pm25_laei2013_breaks)])) #17 (same as breaks)

png('../map_outputs/pm10_2020_sc3_minus_sc1_raster.png', width = 10, height = 8, units = 'in', res = 300)

plot <- levelplot(pm10_2020_sc3_minus_sc1_raster,
          maxpixels = pm10_2020_sc3_minus_sc1_raster@ncols/2 * pm10_2020_sc3_minus_sc1_raster@nrows/2,
          margin = FALSE,
          colorkey = list(
            at = seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = length(pm10_laei2013_breaks)),
            space = 'right',
            labels = list(at=seq(min(pm10_laei2013_breaks), max(pm10_laei2013_breaks), length = length(pm10_laei2013_breaks)), 
                          labels = paste(" \n \n ",pm10_laei2013_labels), 
                          font = 1,
                          cex = 1.5)
          ),
          par.settings = list(
            axis.line =list( col = 'transparent')
          ),
          scales = list(draw = FALSE),
          col.regions = pm10_laei2013_colours,
          at = pm10_laei2013_breaks)

plot

dev.off()

## png 
##   2

print(plot)

rm(pm10_2020_sc3_minus_sc1_raster, plot)