In christiaanpauw/AQoffset: Makes calculation needed for air quality offsetting

knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE)

  wd <- "//Users/christiaanpauw/Dropbox (Nova SA)/Nova House/Business units/Data-DAT/Controlled-DAT/Rpackages/Rpackages/AQoffset/"
  setwd("/Users/christiaanpauw/Dropbox (Nova SA)/Nova House/Business units/Data-DAT/Controlled-DAT/Rpackages/Rpackages/AQoffset/")
  rdir <-  file.path(paste(wd, "R/", sep=""))
  datadir <- normalizePath(paste("Data/", sep = ""))
  graphdir <- normalizePath(paste("Graph/", sep = ""))
  datadir2 <- "~/Documents/AQ_DATA/Data/"

library(AQoffset)
library(AQoffsetData)
names(hh_24) <- gsub("_pm10_24h", "_24h_pm10", names(hh_24))
names(hh_24) <- gsub("_so2_24h", "_24h_so2", names(hh_24))

wgs.proj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
new.proj <- "+proj=utm +south +zone=35 +ellps=WGS84" 
kwa.proj <- "+proj=utm +zone=35 +south +datum=WGS84 +units=km +ellps=WGS84 +towgs84=0,0,0"

kwazaUTM <- spTransform(kwaza, CRS("+proj=utm +zone=35 +south +datum=WGS84 +units=km +ellps=WGS84 +towgs84=0,0,0"))
kwazaUTM$n <- 1
kwazaUTM <- kwazaUTM[,grep("n", names(kwazaUTM@data))]

#kr <- rasterize(kwazaUTM, r3, "n", fun="count")
#levelplot(kr) + latticeExtra::layer(sp::sp.polygons(kwazaUTM, lwd = 0.25))

#remove the other datasets here to save memory
#load(paste(datadir, "/kwaza_eskom_24h_camx.Rda", sep = ""))
#load(paste(datadir, "/eskom.Rda", sep = ""))
#load(paste(datadir, "/kwaza_eskom_year_camx.Rda", sep = ""))
#load(paste(datadir, "/kwaza_household_energyRaster.Rda", sep="")) # households, total households

load(paste(datadir, "/EOP-fuel_users.Rda", sep="")) # dta
kpext <- extent(KWA)
ylabel <- expression(paste("Distribution of daily average ", mu, g/m^{3}, "over 364 days"))
hylabel <- expression(paste("Distribution of daily average ", mu, g/m^{3}, "over 182 days"))
refras <- raster(extent(all_base), ncol = 100, nrow = 100)
eskom_kwa_base <- mask(eskom_kwa_base, masker)
eskom_kwa_base <- resample(eskom_kwa_base, refras)
eskom_kwa_proj <- mask(eskom_kwa_proj, masker)
eskom_kwa_proj <- resample(eskom_kwa_proj, refras)
hh_24 <- resample(hh_24, refras)
all_base2 <- all_base
all_base <- resample(all_base, refras)
masker2 <- resample(masker, refras)

Baseline scenario

Baseline emissions

Household emission sources

KwaZamokuhle is characterised by a high proportion of households who use coal for domestic cooking and heating. The importance of coal is visible in the results on the question on the main energy carrier for heating from the 2011 Census; these results are shown in the Table 1.

kable(kwa_census_df, row.names = FALSE, caption = "Main energy carrier for heating from the 2011 Census")
cat("\n")

#levelplot(households[[-c(7:8)]], pretty = TRUE, par.settings=BuRdTheme, main = "Number of households by energy carrier \nused for heating: Census 2011", scales=list(draw=FALSE)) + addkwa()

It is clear that there are far fewer users of wood and dung than of coal. Because there are very few households who use wood that do not on occasion use coal as well, the responses for coal on the question "Mark ALL the energy carriers that you use for heating?" are used as an estimate for domestic solid fuel use. A summary is provided in Table 2.

kable(kwa_census_fuel_df, row.names = FALSE, caption = "Solid fuel use for heating from the 2011 Census")

It is known however that because the Census asks only a question about the main energy carrier, the number of solid fuel users are underestimated. The results of the survey lead to a substantially higher estimate compared to the census. The estimates of coal using households per sub-place derived from the household survey are shown in Table 3. It is clear that the number of coal using households derived from the survey results is substantially higher.

kable(data.frame(coal.heat[!is.na(coal.heat[,2]), ], row.names = NULL), caption = "Estimated number of coal using households per subplace with upper and lower bound of the 95% confidence interval")

The approximate spatial distribution of the households by coal use are shown below. The high concentration of coal users in the southern and eastern parts of KwaZamokuhle is clear.

sbcoal <- paste("Count (total of ",cellStats(coalPop, sum), ") following Census 2011", sep ="")

coalPop <- resample(coalPop, refras)
coalPop[coalPop == 0] <- NA
levelplot(crop(coalPop, kpext), par = BuRdTheme, margin = NA, at = maak_at(coalPop, beginnul = TRUE), main = "Estimated location of coal using household in Kwazamokuhle ", scales=list(draw=FALSE), sub = sbcoal) + addkwa() + txtkwa()

Baseline emissions are calculated from the results of a domestic fuel use survey. The estimates for fuel consumption based on the household survey are shown in Table 4. Winter coal consumption is understandably higher than summer consumption. Once again KwaZamokuhle SP has the highest consumption.

names(dta) <- gsub("suburb", "", names(dta))
names(dta) <- gsub("households", "#HH", names(dta))
names(dta) <- gsub("month.kg", "kg/m", names(dta))
names(dta) <- gsub("perc", "%", names(dta))
names(dta) <- gsub("winter.ave", "ave(W)", names(dta))
names(dta) <- gsub("summer.ave", "ave(S)", names(dta))
names(dta) <- gsub("\\.", " ", names(dta))
row.names(dta) <- NULL
pander(dta, split.cells = 40, 
       split.table = Inf, 
       rownames = NULL, 
       caption = "Baseline coal use by Suburb")

Baseline states

The modelled baseline PM10 and SO~2~ resulting from household emissions in KwaZamokuhle is shown below. The tiles represent the 25th and 50th percentiles, the mean and the 75th and 99th percentiles. The last tile in each row represent the standard deviation.

hh_24_resample <- resample(hh_24, refras)

raster_dist_plot(hh_24_resample, multi = c("pm10", "so2"), th=BuRdTheme, mn = "Baseline states from households", sb = hylabel) + addkwa() + txtkwa()

The distributions of the PM10 and SO~2~ are fairly similar but the concentration of PM10 is modelled to be higher. It is also clear that the concentration of both PM10 and SO~2~ decreases rapidly with distance from the household sources.

\pagebreak The modelled baseline secondary PM2.5 resulting from the Hendrina power station is shown below.

raster_dist_plot(eskom_kwa_base, multi = c("pm10", "so2"), mn = "Secondary sulphate particulates and SO2 from Hendrina PS\n in the baseline scenario in the region of KwaZamokuhle", sb = ylabel, verbose = TRUE) + txtkwa() + addkwa() 

kable(raster_dist_plot(eskom_kwa_base, multi = c("pm10", "so2"), tab_out = TRUE), digits = 2, caption = "Secondary sulphate particulates and SO2 from Hendrina PS in the baseline scenario in the region of KwaZamokuhle (annual ug/m3)")

Exceedence of specified concentrations in the baseline scenario

The count of days where the PM10 concentrations that resulted from household emissions and the industrial point source is modelled to exceed a specified level, is shown below for the baseline scenario.

hh_pm10_exceed = count_exceed(hh_24, pol = "pm10", min = 50, max = 300, by = 25, knip = FALSE)
rasterVis::levelplot(hh_pm10_exceed, main = "Count of days when PM10 from households \nexceeded specified concentration", sub = "182 days", par.settings=viridisTheme(), at = maak_at(hh_pm10_exceed, beginnul = TRUE)) + addkwa()

The exceedances of the daily PM10 standard related to household emissions occur over the southern and eastern side of KwaZamokuhle with the highest exceedance count at every level occurring in KwaZamokuhle SP.

hh_so2_exceed = count_exceed(hh_24, pol = "so2", min = 50, max = 300, by = 25, knip = FALSE)
rasterVis::levelplot(hh_so2_exceed, main = "Count of days when SO2 from households \nexceeded specified concentration", sub = "182 days", par.settings=viridisTheme(), at = maak_at(hh_so2_exceed, beginnul = TRUE)) + addkwa() + txtkwa()
 ```

As far as the Industrial point source is concerned, the baseline scenario is the compliance scenario. For this reasons the baseline emissions and the associated ambient concentrations are very low. 

```r
eskom_pm10_exceed_base = count_exceed(eskom_kwa_base[[grep("pm", names(eskom_kwa_base))]], pol = NULL, min = 0, max = 20, by = 1, knip = TRUE)

eskom_pm10_exceed_proj = count_exceed(eskom_kwa_proj[[grep("pm", names(eskom_kwa_proj))]], pol = NULL, min = 1, max = 30, by = 5, knip = TRUE)

levelplot(eskom_pm10_exceed_base, main = "Count of days when concentrations of secondary PM from \nEskom exceeded specified concentration in the baseline scenario", sub = "364 days", scales=list(draw=FALSE), par.settings=viridisTheme, at = maak_at(eskom_pm10_exceed_base, beginnul = TRUE)) + addkwa()

count_exceed_table(eskom_pm10_exceed_base, cp = "Count of days when concentrations of secondary PM from \nEskom exceeded specified concentration in the baseline scenario")

The count of days where SO~2~ from the industrial point source exceeded the specified level for the baseline scenario is shown below.

eskom_so2_exceed_base = count_exceed(eskom_kwa_base[[grep("so2", names(eskom_kwa_base))]], pol = NULL, min = 0, max = 10, by = 1, knip = TRUE)

eskom_so2_exceed_proj = count_exceed(eskom_kwa_base[[grep("so2", names(eskom_kwa_proj))]], pol = NULL, min = 0, max = 10, by = 1, knip = TRUE)

levelplot(eskom_so2_exceed_base, main = "Count of days when secondary SO2 from \nEskom exceeded specified concentration in the baseline scenario", sub = "364 days", scales=list(draw=FALSE), par.settings=viridisTheme, at = maak_at(eskom_so2_exceed_base, beginnul = TRUE)) + addkwa() + txtkwa()

#bar_exceed(eskom_pm10_exceed, ttl = "Aggregated of days when PM10 from Eskom \nexceeded specified concentration")
count_exceed_table(eskom_pm10_exceed_base, cp = "Baseline exceedence count: PM form Eskom")

\pagebreak The count of days where SO~2~ concentration that resulted from household emissions is modelled to exceed a specified level, is shown below for the baseline scenario.

hh_so2_exceed = count_exceed(hh_24, pol = "so2", min = 100, max = 300, by = 25, knip = FALSE)

levelplot(hh_so2_exceed, main = "Count of days when SO2 from households \nexceeded specified concentration", sub = "182 days", par = viridisTheme(), scales=list(draw=FALSE), at = maak_at(hh_so2_exceed, beginnul = TRUE)) + addkwa() + txtkwa()

Occurrences of exceedances of the daily SO~2~ standard related to coal combustion from households are less common than that for PM10. The area over which such occurrences are smaller than that of PM10 for every level and the count of exceedances for each level is also less. Once again most exceedances occur in KwaZamokuhle SP.

#bar_exceed(hh_so2_exceed, ttl = "Aggregated count of days when SO2 from \nhouseholds exceeded specified concentration")
count_exceed_table(hh_so2_exceed)

Combined baseline states

# create a people raster of the same extent as hh_24
households <- rasteriseCensus(KWA, ref = all_base, refres = c(nrow(all_base), ncol(all_base)), verbose = TRUE)
total_households <- calc(households, sum, na.rm = TRUE)
names(total_households) <- "Total_households"
total_households[total_households == 0] <- NA
people <- total_households *3.7
people[people==0] <- NA
names(people) <- "total"
all.equal(extent(people) , extent(all_base))

raster_dist_plot(mask(all_base, masker2), multi = c("pm10", "so2"), mn = "Mean 24h baseline concentrations \nin the baseline scenario from households and Hendrina PS", sb = hylabel) + addkwa() + txtkwa()

raster_dist_plot(all_base, multi = c("pm10", "so2"), meanonly = TRUE, mn = "Mean 24h baseline concentrations \nin the baseline scenario from households and Eskom", sb = hylabel) + addkwa() + txtkwa()

kable(raster_dist_plot(all_base, multi = c("pm10", "so2"),tab_out  = TRUE), caption = "Mean 24h baseline concentrations \nin the baseline scenario from households and Eskom (ug/m3 over 182 days)", digits = 2)

Project boundary

The spatial extent of the project boundary is the overlapping extent of the ambient contribution of the baseline emissions above 2 ug/m3 per year or 19 ug/m3 per day in PM10 or SO~2~, and the same for project emissions from the managed activity. Project emissions from the managed activity are the emissions from Eskom in the business-as-usual scenario. Baseline emissions from households are the emissions from households before implementation of the intervention.

Application of thresholds to determine project boundary

The baseline emissions from household coal use are summarised below. The tiles represent the 1st, 25th and 50th percentiles, the mean and the 75th and 99th percentiles. The last three tiles represent the standard deviation, inter-quartile range and 99% range of the data.

raster_dist_plot(all_base2, multi = c("pm10", "so2"), mn = "Definition of the project boundary", sb = hylabel, th=BuRdTheme, verbose = TRUE) + addkwa() + txtkwa()

\pagebreak The extent of the impact of household emissions in the baseline scenario above an annual average of 2ug/m3 or an daily average of 19ug/m3 is shown below. Blank cells are outside of the project boundary.

all_base <- resample(project_boundary(all_base2), refras)
raster_dist_plot(all_base, multi = c("pm10", "so2"), meanonly = FALSE, mn = expression(paste("Definition of the project boundary: 24h ", mu,g/m^{3})),  sb = "Non-blank areas are where daily or \n annual PM10 exceed the threshold (182 days)") + addkwa() + txtkwa()

Summary

raster_dist_plot(all_base, multi = c("pm10", "so2"), meanonly = TRUE, mn = expression(paste("Definition of the project boundary: 24h ", mu,g/m^{3})),  sb = "Non-blank areas are where daily or \n annual PM10 exceed the threshold (182 days)", verbose = TRUE) + addkwa() + txtkwa()

It is clear from the application of the threshold that, if the modelling is correct, the impact of domestic burning is localised in close proximity to the emissions. The project boundary includes the whole of the main place KwaZamokuhle. Of all the sub-places that make up the main place KwaZamokuhle, the highest mean and maximum concentrations are found in KwaZamokuhle SP.

Baseline impact

Baseline impact can be determined using different calculation approaches. Four approaches will be demonstrated. These are: 1. Health risk approach 2. Particle equivalence approach 3. Standards weighted intake 4. Burden of disease approach

Health risk approach

The health risk approach uses an air quality index based on the relative risk of short term mortality associated with every pollutant. Here the pollutants are PM10 and SO~2~, but O~3~ and NO~2~ can potentially be added. The baseline impact represents the combined impact of all pollutants. A summary of the daily API resulting from households is given below.

kwaza_API <- rasterAPI(all_base, idpos = 1, polpos = 2, cyclepos = 5)

raster_dist_plot(kwaza_API,  multi = NULL,  mn = "Distribution of baseline API \nfor the baseline scenario", sb = "Air pollution index over 182 days") + addkwa() + txtkwa()


raster_dist_plot(kwaza_API,  multi = NULL,  mn = "Distribution of baseline API \nfor the baseline scenario", sb = "Air pollution index over 182 days", tab_out = TRUE) 
kable(raster_dist_plot(kwaza_API, tab_out = TRUE, multi = NULL), caption = "Distribution (mean of percentiles) of baseline API \nfor the baseline scenario", digits = 2, col.names = "AP Index value")

kable(raster_dist_plot(kwaza_API, meanonly = TRUE,
                 multi = NULL, 
                 mn = "Distribution of baseline API \nfor the baseline scenario",
                 sb = "Mean index value over 182 days", 
                 tab_out = TRUE), digits = 2, col.names = "AP Index value")

\pagebreak The count of days where the baseline API is modelled to exceed a specified level is shown below for the baseline scenario.

levelplot(count_exceed(kwaza_API, min = 3, max = 10, by = 1, knip = FALSE) , 
          par = viridisTheme, 
          scales=list(draw=FALSE), 
          main = "Days in KwaZamokuhle where baseline API exceeded specified level", 
          sub = "182 days") + addkwa() + txtkwa()

KwaZamokuhle SP has the highest health risk related to short term exposure.

Particle equivalence approach

With sources that have a very localised dispersion, the particle equivalent impact (as PM10 equivalent) of that source is simply the concentration of the PM10 emitted by that source at every receptor; therefore, in the case of household emissions, the baseline impact is the same as the baseline state of PM10 shown above.

The contribution of households and Eskom to particulate concentrations are shown below, first in detail and thereafter only the contribution to the 182 day mean.

```{R combined_base} raster_dist_plot(combined_base, multi = c("Households_pm10_24h\.","Eskom_24h_pm10\.", "All_pm10_24h\."), meanonly = FALSE, mn = "Baseline PM contribution", sb = hylabel) + addkwa() + txtkwa()

### Standards weighted intake

The standards weighted intake approach is roughly equivalent to a population weighted air pollution index with the index weights determined by the standard. In this way it is related to the health risk approach that also requires population data. 

```r
people <- crop(people, refras)
people <- resample(people, refras)
swi_base <- SWI(all_base, idpos = 1, polpos = 2, aveperiodpos = 3, cyclepos = 5, pop = people)


ppatt <- seq(0, ceiling(maxValue(people)), ceiling(maxValue(people))/20)
#1 Pop density
p1 <- rasterVis::levelplot(people, par.settings = BuRdTheme, main = "Population density", margin = FALSE, scales = list(draw = FALSE), at = ppatt) + addkwa() + txtkwa()

#2 SWI base
p2att <- seq(0, ceiling(maxValue(swi_base)), ceiling(maxValue(swi_base))/20)
p2base <- levelplot(swi_base, par.settings = BuRdTheme, scales = list(draw = FALSE), main = "SWI: baseline scenario", margin = FALSE, at = p2att) + addkwa() + txtkwa()

som <- function(x) {sum(x, na.rm = TRUE)}
all_base_so2_sum <- calc(all_base[[grep("so2", names(all_base))]], som)

#3 CUM SO2
att1 <- seq(0, ceiling(maxValue(all_base_so2_sum)), ceiling(maxValue(all_base_so2_sum))/20)
all_base_so2_sum <- crop(all_base_so2_sum, all_base)
p3base <- levelplot(all_base_so2_sum, par.settings = BuRdTheme, main = expression(paste("Cumilative baseline ", SO[2])), margin = FALSE, scales = list(draw = FALSE), at  = att1) + addkwa() + txtkwa()

all_base_pm_sum <- calc(all_base[[grep("pm10", names(all_base))]], som)

#4 CUM PM
p4att <- seq(0, ceiling(maxValue(all_base_pm_sum)), ceiling(maxValue(all_base_pm_sum))/20)
p4base <- levelplot(all_base_pm_sum, par.settings = BuRdTheme, main = expression(paste("Cumilative baseline ", PM[10])), margin = FALSE, scales = list(draw = FALSE), at = p4att) + addkwa() + txtkwa()

grid.arrange(p3base, p4base, p1, p2base,  nrow=2)

Burden of disease approach

The burden of disease approach quantifies the actual or expected incidence of adverse health outcomes attributible to the exposure to ambient air pollution, and expresses the impact of the air pollution in terms of the proportion or number of cases of a specific outcome or as a weighted aggregate of such outcomes.

The method requires estimates of the incidence rates as well as the total susseptible population of each of the outcomes to be quantified. The method furthermore requires concentration response functions for every outcome. If a monatary estimate of the impact is to be made, a cost estimate per case is also needed. In some cases the concentration response fucntion is defined for PM2.5 but exposure data is availible for PM10. I such cases one needs to make an assumption of the fraction of PM10 that is PM2.5. Chronic and acute effects have to be treated seperately.

kable(summarise.sicklist(endlist.eu)[,c(1:5, 12:15)],digits = 2, caption = "Assumptions for burden of disease calculations")

\pagebreak The baseline burden of disease impact for PM10 in the baseline scenario is shown below.

# estimate incidences
rBoDPM10 <- rasterCREP(sl = endlist.eu,
                       pollutant = "PM10", 
                       cconc = all_base[[grep(pattern = "pm10", x = tolower(names(all_base)),
                                              fixed = TRUE)]], 
                       ppopr = people, 
                       bbase.conc = 10, 
                       rrisk.only = FALSE, 
                       verbose = FALSE, 
                       output = "AM",
                       ddelta = 1)

rBoDPM10 <- aggregateCREPstack(rBoDPM10)
# tabulate total incidences per health outcome
kable(sumIncidences(rBoDPM10), digits = 2)

# plot incidences
levelplot(crop(rBoDPM10, kpext), scales = list(draw = FALSE),
main = "Burden of disease for PM10 in the baseline scenario,
given as number of cases per health outcome", 
          sub = "182 days",
          par = BuRdTheme) + addkwa() + txtkwa()

\pagebreak The baseline burden of disease impact for SO~2~ resulting from household emissions is calculated using concentration-response functions derived from studies in China.

kable(summarise.sicklist(endlist.china)[,c(1:5, 12:15)],digits = 2, caption = "Assumptions for burden of disease calculations")

Results are shown below

rBoDSO2 <- rasterCREP(sl = endlist.china,
                      pollutant = "SO2", 
                      cconc = all_base[[grep(pattern = "so2", 
                                           x = tolower(names(all_base)), 
                                           fixed = TRUE)]], 
                      ppopr = people, 
                      bbase.conc = 10, 
                      rrisk.only = FALSE,
                      output = "AM",
                      verbose = FALSE,
                      ddelta = 1)

# summarise the number of cases across all days
rBoDSO2 <- aggregateCREPstack(rBoDSO2)

# tabulate total incidences per health outcome
kable(sumIncidences(rBoDSO2), digits = 2)

levelplot(crop(rBoDSO2, kpext),
          main = "Burden of disease for SO2 
from all sources,
given as # of cases per health outcome",
          sub = "182 days",
          par = BuRdTheme) + addkwa() + txtkwa()

\pagebreak

Project scenario

The project scenario is the implementation of a stove exchange for all RDP houses who use coal, where the households exchange their old coal stoves for a full retrofit and LPG.

The estimates of fuel users per house type are shown below.

kable(fuel_house[!is.na(fuel_house[,2]), ])
fh = fuel_house[which(fuel_house[,2] == "Fuel_RDP"), ]
fh$place = unique(fuel_house[,1])[c(1,3:6)]

The target for the project activity is therefore between r colSums(fh[,3:5])["Lower"] and r colSums(fh[,3:5])["Upper"] with a point estimate of r colSums(fh[,3:5])["PointEst"]

kable(fh, caption = "Implementation targets per sub-place")

Project emissions

Improvement per household

Switch to LPG represents, for all prcticalpurposes, a 100% reduction in PM10 emission.

Kitchen King = 0.018 * Union * 1.4 * 2 to be conservative ~ 5%

Proportion of solid fuel using households reachable by the project

Project states

Project states from households

hh_24_proj <- resample(hh_24, refras) * 0.48
raster_dist_plot(hh_24_proj, mn = "Project states from households ", sb = "182 days", multi = c("pm10", "so2"), th=BuRdTheme) + addkwa() + txtkwa()

The number of high concentration events resulting from household emissions are modelled to decreased in the project scenario.

hh_24_proj_exceed_pm <- count_exceed(hh_24_proj, pol = "pm", min = 50, max= 150, by =25)

levelplot(hh_24_proj_exceed_pm, main = "Count of days when PM from households \nexceeded specified concentration in the project scenario", sub = "182 days", par = viridisTheme(), scales=list(draw=FALSE), at = maak_at(hh_24_proj_exceed_pm, beginnul = TRUE)) + addkwa() + txtkwa()

count_exceed_table(hh_24_proj_exceed_pm)

Project states from Eskom

A summary of the atmospheric states resulting from Eskom emissions in the project scenario over a period of 364 days is shown below.

raster_dist_plot(eskom_kwa_proj, mn = "Modelled project states from Hendrina PS", sb = ylabel, multi = c("pm10", "so2"), th=BuRdTheme) + addkwa() + txtkwa()

Eskom is modelled to contribute very little to high PM concentrations with most days modelled to have less than 6ug/m3

count_exceed_table(eskom_pm10_exceed_proj, cp = "Project exceedance count: PM from Eskom. Daily standard is 75")

SO~2~ from Eskom is also very low with no day modelled to have more than 5 ug/m3 of SO~2~ as a result of Eskom emissions

count_exceed_table(eskom_so2_exceed_proj, cp = "Project exceedance count: SO~2~ from Eskom. Daily standard is 125")

Project states (all sources)

In the project scenario, PM and SO2 emissions from households would be approximately 49% of the baseline scenario. Eskom emissions of SO~2~ in the project scenario is the business-as-usual emissions which is approximately 6.7 times higher than emissions in the compliance scenario (which is the baseline scenario). Since household emissions were only availible for 182 days, the combined comparison is only done for 182 days.

raster_dist_plot(hh_24_proj, mn = "Modelled project states from Households", sb = ylabel, multi = c("pm10", "so2"), th=BuRdTheme) + addkwa() + txtkwa()

allPM10_proj <- stack(hh_24_proj[[grep("pm10", names(hh_24_proj))]], eskom_kwa_proj[[grep("pm10", names(eskom_kwa_proj))]])
yd <- gsub("Eskom_24h_pm10\\.|Households_24h_pm10\\.", "", names(allPM10_proj))
corrdays <- yd[which(duplicated(yd))]
allPM10_proj <- stackApply(allPM10_proj, indices = match(yd , corrdays), sum)
names(allPM10_proj) <- gsub("index_", "All_pm10_24h.", names(allPM10_proj))

mn.pm.proj = expression(paste("Daily average ",PM[10], " in the project scenario"))
sb.pm.proj = expression(paste("24h average ", mu, g/m^{3}, " "))
raster_dist_plot(ss = stack(hh_24_proj, eskom_kwa_proj, allPM10_proj), multi = c("Households", "Eskom", "All"), mn = mn.pm.proj, sb = sb.pm.proj) + addkwa() + txtkwa()

allSO2_proj <- stack(hh_24_proj[[grep("so2", names(hh_24))]], eskom_kwa_proj[[grep("so2", names(eskom_kwa_proj))]])
yd.so2 <- gsub("Eskom_24h_so2\\.|Households_24h_so2\\.", "", names(allSO2_proj))
corrdays.so2 <- yd.so2[which(duplicated(yd.so2))]
allSO2_proj <- stackApply(allSO2_proj, indices = match(yd.so2 , corrdays.so2), sum)
names(allSO2_proj) <- gsub("index_", "All_so2_24h.", names(allSO2_proj))

The combined SO~2 for households and Hendrina Power station for the first r{nlayers(allSO2_proj)} days are shown below with the combined baseline states in row one and the SO~2 from Hendrina Power Station and households in rows two and three respectively.

mn.so2.proj = expression(paste("Daily average ",SO[2], " in the project scenario"))
sb.so2.proj = expression(paste("24h average ", mu, g/m^{3}, " "))
raster_dist_plot(ss = stack(allSO2_proj, hh_24_proj[[grep("so2", names(hh_24))]], eskom_kwa_proj[[grep("so2", names(eskom_kwa_proj))]]), multi = c("Households", "Eskom", "All"), mn = mn.so2.proj, sb = sb.so2.proj) + addkwa() + txtkwa()

all_proj <- stack(allPM10_proj, allSO2_proj)
didx <- grep("NA", names(all_proj))
if (length(didx) > 0) all_proj <- all_proj[[-didx]]
combined_proj <- stack(hh_24_proj, eskom_kwa_proj, all_proj)

mn.all.proj <- expression(paste("Daily average ", PM[10], " and ", SO[2], " \nin the project scenario (all sources)"))

raster_dist_plot(all_proj, multi = c("pm10", "so2"), beginnul = TRUE, mn = mn.all.proj, sb = hylabel) + addkwa() + txtkwa()

To summarise the project states, the mean concentrations resulting from households and Hendrina Power Station as well as the mean combined project state is shown below.

Project impacts

Project standards weighted intake

swi_proj <- SWI(all_proj, idpos = 1, polpos = 2, aveperiodpos = 3, cyclepos = 5, pop = people)

swi_proj_p <- resample(swi_proj, refras)
pp <- resample(people, refras)
swi_proj_stack <- stack(pp , swi_proj_p )
names(swi_proj_stack ) <- c("People", "Std.Weighted Intake")

att1 <- seq(0, ceiling(maxValue(pp)), ceiling(maxValue(pp))/20)
p1 <- levelplot(pp, par.settings = BuRdTheme, main = "Population density", margin = FALSE, scales = list(draw = FALSE), at = att1)+ addkwa() + txtkwa()

att2 <- seq(0, ceiling(maxValue(swi_proj_p)), ceiling(maxValue(swi_proj_p))/20)

p2 <- levelplot(swi_proj_p, par.settings = BuRdTheme, main = "Project SWI", margin = FALSE, scales = list(draw = FALSE), at = att2) + addkwa() + txtkwa()

all_proj_so2_sum <- resample(calc(all_proj[[grep("so2", names(all_proj))]], som), y = refras)

p3 <- levelplot(all_proj_so2_sum, par.settings = BuRdTheme, main = expression(paste("Cumilative project ", SO[2])), margin = FALSE, scales = list(draw = FALSE), at = maak_at(all_proj_so2_sum)) + addkwa() + txtkwa()

all_proj_pm_sum <- resample(calc(all_proj[[grep("pm10", names(all_proj))]], som), y = refras)

p4 <- levelplot(all_proj_pm_sum, par.settings = BuRdTheme, main = expression(paste("Cumilative project ", PM[10])), margin = FALSE, scales = list(draw = FALSE)) + addkwa() + txtkwa()

grid.arrange(p4, p3, p1, p2, ncol=2)

Projecth health risk approach

Project impacts

kwaza_API_proj <- rasterAPI(all_proj, idpos = 1, polpos = 2, aveperiodpos = 3, cyclepos = 5)

raster_dist_plot(resample(kwaza_API_proj, y = refras), multi = NULL, mn = "Distribution of project API \nfor the project scenario", sb = "182 days", diverge = TRUE) + addkwa() + txtkwa()

Project burden of disease

rBoDPM10_proj <- rasterCREP(sl = endlist.eu,
                       pollutant = "PM10", 
                       cconc = all_proj[[grep("pm10", names(all_proj))]], 
                       ppopr = people, 
                       bbase.conc = 10, 
                       rrisk.only = FALSE, 
                       verbose = FALSE, 
                       output = "AM")

rBoDPM10_proj <- aggregateCREPstack(rBoDPM10_proj)

levelplot(resample(rBoDPM10_proj, refras), scales = list(draw = FALSE),
          main = "Burden of disease for PM10 in the project scenario,
given as number of cases per health outcome",
sub = "182 days",
par = BuRdTheme) + addkwa() + txtkwa()

# tabulate total incidences per health outcome
kable(sumIncidences(rBoDPM10_proj), digits = 2, caption = "Estimated burden of disease in the project scenario")

rBoDSO2_proj <- rasterCREP(sl = endlist.china,
                       pollutant = "SO2", 
                       cconc = all_proj[[grep("so2", names(all_proj))]], 
                       ppopr = people, 
                       bbase.conc = 10, 
                       rrisk.only = FALSE, 
                       verbose = FALSE, 
                       output = "AM")

rBoDSO2_proj <- aggregateCREPstack(rBoDSO2_proj)

levelplot(resample(rBoDSO2_proj, refras), scales = list(draw = FALSE),
          main = expression(paste("Burden of disease for ", SO[2] ," in the project scenario, given as number of cases per health outcome")),
sub = "182 days",
par = BuRdTheme) + addkwa() + txtkwa()

Impact reduction

Particle precursor approach

In the particle precursor approach, the impact is defined as the difference in particulate concentrations between the baseline and project scenarios

The impact of the project scenario in terms of PM10 concentrations are shown below

PM_impact <- resample(all_base[[grep("pm10", names(all_base))]] - all_proj[[grep("pm10", names(all_proj))]], refras)

raster_dist_plot(PM_impact, mn = "Impact of the project activtiy on PM10 concentrations", sb = "Daily differenence between baseline and project scenario over 182 days", diverge = TRUE) +addkwa() + txtkwa() 

In the same wat as above, the impact of the project scenario on SO~2~ concentrations are shown below.

SO2_impact <- resample(all_base[[grep("so2", names(all_base))]] - all_proj[[grep("so2", names(all_proj))[1:182]]], refras)

raster_dist_plot(SO2_impact, mn = expression(paste("Impact of the project activtiy on ", SO[2], " concentrations (", mu, g/m^{3},")")), sb = "Daily differenence between baseline and project scenario over 182 days", diverge = TRUE) +addkwa() + txtkwa() 

The combined impact of the project activity on PM~10 and SO~2 emissions are shown below

combined_impact <- stack(PM_impact, SO2_impact)
raster_dist_plot(combined_impact, mn = expression(paste("Impact of the project activtiy on ", SO[2], " and ", PM[10], " concentrations (", mu, g/m^{3},")")), multi = c("pm10", "so2"),sb = "Daily differenence between baseline and project scenario over 182 days", diverge = TRUE) +addkwa() + txtkwa()

Health risk approach

api_impact <- kwaza_API - kwaza_API_proj
api_impact <- resample(mask(api_impact, masker2), refras)

raster_dist_plot(api_impact, multi = NULL,  mn = "Distribution of daily API improvement \ndue to the project scenario",   sb = "182 days", diverge = TRUE) + addkwa() + txtkwa()

raster_dist_plot(api_impact, meanonly = TRUE, multi = NULL,  mn = "Mean daily API improvement \ndue to the project scenario",   sb = "182 days", diverge = TRUE) + addkwa() + txtkwa()

Standard weighted intake approach

swi_base[is.na(swi_base)] <- 0
swi_proj[is.na(swi_proj)] <- 0
swi_impact <- resample(swi_base - swi_proj, refras)

swi_impact[swi_impact==0] <- NA
diverge0(levelplot(swi_impact, main = "Impact in terms of reduction \nin standards weighted intake", sub = "Difference between baseline and \nproject scenarion over 182 days", margin = FALSE, par.settings = BuRdTheme), ramp = "RdBu") + addkwa() + txtkwa()

Burden of disease approach

rBoDPM10[is.na(rBoDPM10)]  <- 0
rBoDPM10_proj[is.na(rBoDPM10_proj)] <- 0
rBoDPM10_impact <- rBoDPM10 - rBoDPM10_proj
rBoDPM10_impact[rBoDPM10_impact==0] <- NA

att = seq(0, ceiling(max(maxValue(rBoDPM10_impact))), ceiling(max(maxValue(rBoDPM10_impact)))/ 20)

levelplot(rBoDPM10_impact, 
          main = "Impact in terms of disease burden due to the \n implementation of the project activity", 
          sub = "number of cases", 
          par = BuRdTheme, 
          scales = list(draw = FALSE), at =att) + addkwa() + txtkwa()

kable(sumIncidences(rBoDPM10_impact), digits = 2, caption = "Estimated burden of disease improvement within the project boundary")

Impact (number to implement)

Impact deficit: particle percursor approach

names(eskom_kwa_base) <- paste("Base", names(eskom_kwa_base))
names(eskom_kwa_proj) <- paste("Proj", names(eskom_kwa_proj))

eskomstack <- stack(eskom_kwa_base, eskom_kwa_proj)
names(eskomstack) <- gsub("\\.Eskom_24h","", names(eskomstack))

eskom_deficit <- eskom_kwa_base - eskom_kwa_proj

plabel <- expression(paste("Daily average ", mu, g/m^{3}, " "))

eskom_deficit <- setValues(x = eskom_deficit, values = abs(getValues(eskom_deficit)))

kable(make_impact_table(base = eskom_kwa_base, proj = eskom_kwa_proj, def = eskom_deficit, pmvar = "pm"), caption = paste("Summary of Industrial emissions baseline, project and deficit (in", expression(paste("Daily average ", mu, g/m^{3}, " PM10"))))

kable(make_impact_table(base = eskom_kwa_base, proj = eskom_kwa_proj, def = eskom_deficit, pmvar = "so2"), caption = paste("Summary of Industrial emissions baseline, project and deficit (in", expression(paste("Daily average ", mu, g/m^{3}, " SO2"))))

eskom_deficit[eskom_deficit == 0] <- 0.000001
names(eskom_deficit) <- paste("Deficit", names(eskom_deficit), sep="")
raster_dist_plot(eskom_deficit, multi = c("pm10", "so2"), mn = "Distribution of daily impact deficit", sb = plabel, beginnul = FALSE) +addkwa() + txtkwa()

hh_prop_pp <- eskom_deficit[[grep("pm10", names(eskom_deficit))]] / hh_24[[grep("pm10", names(hh_24))]]

raster_dist_plot(hh_prop_pp, mn = "Proportion of household emission\n required to offset Eskom emissions", maksimum = 0.5, beginnul = TRUE, meanonly = FALSE) + addkwa() + txtkwa()

Impact deficit: health risk approach

The health risk approach

api.deficit <- rasterAPI(eskom_deficit, idpos = 1, polpos = 3, aveperiodpos = 2, cyclepos = 5, per_source = TRUE)

api.prop <- api.deficit  / kwaza_API
if (any(is.na(getValues(api.prop)))) api.prop[is.na(api.prop)] <- 0

raster_dist_plot(api.prop, mn = "Required proporion of household API", sb = "API index points", maksimum = 10) + addkwa() + txtkwa()
mean(cellStats(api.deficit, mean))

Impact deficit: standards weighted mean

raster_dist_plot(eskom_deficit, meanonly = TRUE) + addkwa() + txtkwa()

swi.deficit <- SWI(eskom_deficit)
swi.hh <- SWI(hh_24)
swi.prop <- swi.deficit  / swi.hh
swi.deficit <- resample(swi.deficit, refras)

diverge0(levelplot(crop(swi.deficit, kpext), main = "Impact deficit in the SWI approach", sub = "Standards weighted intake", scales = list(draw = FALSE), margin = FALSE, par = BuRdTheme, at = maak_at(swi.deficit, beginnul = TRUE)), ramp = "RdBu") + addkwa() + txtkwa()

levelplot(crop(resample(swi.prop, refras),kpext), main = "Required proporion of household SWI", sub = "Standards weighted intake", margin = FALSE, par = BuRdTheme, scales = list(draw = FALSE)) + addkwa() + txtkwa()

mean(cellStats(swi.deficit, sum))

Impact deficit: BoD

bod_deficit <- rasterCREP(sl = endlist.eu,
                       pollutant = c("pm10","so2"), 
                       cconc = eskom_deficit, 
                       ppopr = people, 
                       bbase.conc = 0, 
                       rrisk.only = FALSE, 
                       verbose = FALSE, 
                       output = "AM")

bod_deficit <- aggregateCREPstack(bod_deficit)

att = seq(0, max(maxValue(bod_deficit)),  max(maxValue(bod_deficit))/20)
levelplot(resample(bod_deficit, refras), main = "Burden of disease deficit",
sub = "Cases over 182 days. Reference concentration: 0",
par = BuRdTheme, at = att) + addkwa() + txtkwa()

kable(sumIncidences(bod_deficit), digits = 2)

bod.hh <- rasterCREP(sl = endlist.eu,
                       pollutant = c("pm10", "so2"), 
                       cconc = resample(crop(hh_24, eskom_deficit), eskom_deficit),
                       ppopr = people, 
                       bbase.conc = 10, 
                       rrisk.only = FALSE, 
                       verbose = FALSE, 
                       output = "AM")

bod.hh <- aggregateCREPstack(bod.hh)

#pdf(file = "Graph/bod.hh.p.pdf", width = 14)
bod.hh.p <- rasterVis::levelplot(bod.hh, main = "Burden of disease households", sub = "Cases over 182 days. Reference concentration: 10",  par = RdBuTheme, at = maak_at(bod.hh, simetries = TRUE)) + addkwa() + txtkwa()
bod.hh.p
#dev.off()

kable(sumIncidences(bod.hh), digits = 2)

# bod.prop <- bod_deficit  / bod.hh
# 
# maks <- function(x) max(x, na.rm = TRUE)
# 
# bod.prop <- calc(bod.prop, maks)
# 
# levelplot(d.prop, main = "Required % of household required \nto offset the impact of the managed activity (BoD apporach)", sub = "Disease burden", margin = FALSE, par = BuRdTheme) + addkwa() + txtkwa()
# 
# kable(sumIncidences(bod.prop), digits = 2, caption = "BoD improvement within the project boundary due to the project activity")

christiaanpauw/AQoffset documentation built on May 13, 2019, 6:57 p.m.