Skip to content
Snippets Groups Projects
Commit c90866b1 authored by Ben Anderson's avatar Ben Anderson
Browse files

moved to SAVE

parent 6832dd6e
No related branches found
No related tags found
No related merge requests found
Showing
with 0 additions and 9307 deletions
# test SAVE SDRC 2.1 baseline IPF (based on time-use) against
# SDRC 2.2 based on observd kWh
# Libraries ----
library(data.table)
library(drake)
library(ggplot2)
library(here)
library(hms)
library(lubridate)
library(survey)
# Paramaters ----
pPath <- paste0(here::here(), "/plots/") # where to put plots
cPath <- "~/Data/UK_Census/Census_LUTS/2011/EngWales/"
# Rmd
rmdRoot <- paste0(here::here(), "/analysis/compareModels/")
# Load Census data (LSOAs) ----
# LSOA 2001 <-> LSOA 2011 LUT
lsoaLutDT <- data.table::fread(paste0(cPath,
"LSOA01_LSOA11_EW_LU.zip Folder/LSOA01_LSOA11_EW_LUv2.csv"))
sotonLSOA_01DT <- lsoaLutDT[LSOA01NM %like% "Southampton"]
sotonLSOA_11DT <- lsoaLutDT[LSOA11NM %like% "Southampton"]
# unchanged?
table(sotonLSOA_01DT$CHGIND)
table(sotonLSOA_11DT$CHGIND)
sotonLSOA_01DT[, zonecode := LSOA01CD] # set to 2001 LSOA for merge
setkey(sotonLSOA_01DT, zonecode)
sotonLSOA_11DT[, zonecode := LSOA11CD] # set to 2011 LSOA for merge
setkey(sotonLSOA_11DT, zonecode)
# Load Census data (OAs) ----
# OA 2001 <-> OA 2011 LUT
oa2001_2011LutDT <- data.table::fread(paste0(cPath,
"OA01_OA11_LAD11_EW_LU_csv/OA01_OA11_LAD11_EW_LU.csv"))
# not used - likange done at LSOA level
# Load Cesus data (OA <-> LSOA) ----
oa2011LutDT <- data.table::fread(paste0(cPath,
"OA11_LSOA11_MSOA11_LAD11_EW_LU.zip Folder/",
"OA11_LSOA11_MSOA11_LAD11_EW_LUv2.csv"))
sotonOA_11DT <- oa2011LutDT[LAD11NM %like% "Southampton"]
# Function: SDRC 2.1 Time Use model ----
getSDRC2.1 <- function(version){
message("Loading: ", version)
sdrc2.1tuPath <- "~/University of Southampton/SAVE - Documents/WP2-Customer-Model/SDRC-v2.1/results/ONS-TU-2000/"
sdrc2.1f <- paste0("ONS-TU-2000-SAVE-SDRC-v2.1-household-diary-file-weekdays-winter-",version,".dta")
# get the data
message("Loading: ", sdrc2.1f)
df <- paste0(sdrc2.1tuPath,sdrc2.1f)
dt <- data.table::as.data.table(foreign::read.dta(df))
# convert to kWh
dt[, power_sumkW := power_sum/1000]
dt[, power_sumkWh := power_sumkW/2] # cos it's half hours
setkey(dt, hhid_numeric)
uniqueN(dt$hhid_numeric)
# add weights
# weights files
sdrc2.1wPath <- "~/University of Southampton/SAVE - Documents/WP2-Customer-Model/SDRC-v2.1/results/ONS-TU-2000-sms/data/"
#sdrc2.1wf <- "output/ONS-TU-2000-Census-2001-wtdh_ug-Southampton-simoutput-weights.dta"
sdrc2.1wf <- "output/ONS-TU-2000-Census-2001-wtdh_ug-Southampton-v3m1_weekdays_winter_simoutput_weights.dta"
message("Loading weights: ", sdrc2.1wf)
df <- paste0(sdrc2.1wPath,sdrc2.1wf)
wdt <- data.table::as.data.table(foreign::read.dta(df))
setkey(wdt, hhid_numeric)
uniqueN(wdt$hhid_numeric)
# join them
jdt <- wdt[dt[, .(s_halfhour, ba_hrpindsurvey, hhid_numeric, power_sum, power_sumkWh)], allow.cartesian=TRUE]
summary(head(jdt))
# aggregate to LSOA (zonecode = LSOA)
# set survey design so we can do weighted stats
message("Setting svydesign ", version)
svy.jdt <- survey::svydesign(ids = ~hhid_numeric,
weights = ~weight,
data = jdt[!is.na(weight)])
message("Aggregating ", version)
svyt <- survey::svyby(~power_sumkWh, ~s_halfhour+zonecode, svy.jdt, # this is where you'd put the survey categories to average over
svymean,
keep.var = TRUE,
exclude=NULL,
na.action=na.pass
)
dt <- as.data.table(svyt)
dt[, dateTime := lubridate::as_datetime(s_halfhour, tz = "UTC")]
dt[, hms := hms::as.hms(dateTime)]
dt[, mean_powerkWh := power_sumkWh]
setkey(dt, zonecode)
setkey(sotonLSOA_01DT, zonecode)
# > plot LSOAs ----
plotDT <- dt[sotonLSOA_01DT, allow.cartesian=TRUE] # just Soton
p <- ggplot2::ggplot(plotDT,
aes(x = hms,
y = mean_powerkWh,
group = zonecode,
colour = zonecode)) +
geom_line()
p <- p + theme(legend.position = "None") + # no legend
labs(x = "Time of Day",
y = "Mean household kWh per half hour",
caption = paste0("2001 Census & Time Use Model (Model ", version," ): Winter weekdays\n",
uniqueN(plotDT$zonecode), " Southampton LSOAs")
)
ggplot2::ggsave(paste0(pPath,"lsoaProfiles_Winter-Weekdays-Time-Use-Model-", version, "-Sim.png"),
p)
return(dt)
}
# Function: SDRC 2.2 Event day ----
# Can load from OA file
# useful for spurious outier effects discussion
getSDRC2.2ED <- function(){
sdrc2.2Path <- "~/Data/SAVE/processed/"
f <- "SAVE_SDRC_2.2_OA_Event_Day_kWh_profiles_v1.csv.gz"
sdrc2.2_ED_OA_dt <- data.table::fread(paste0(sdrc2.2Path,
f))
# collapse to LSOAs
sdrc2.2_ED_OA_dt[, zonecode := oaCode]
setkey(sdrc2.2_ED_OA_dt, zonecode)
sotonOA_11DT[, zonecode := OA11CD]
setkey(sotonOA_11DT, zonecode)
sdrc2.2_ED_OA_dt <- sotonOA_11DT[sdrc2.2_ED_OA_dt]
dt <- sdrc2.2_ED_OA_dt[, .(meankWh = mean(meanHalfHourkWh)),
keyby = .(LSOA11CD, obsDateTime30m, obsDate)]
dt[, obsDateTime := lubridate::as_datetime(obsDateTime30m)]
dt[, hms := hms::as.hms(obsDateTime, tz = "UTC")] # just in case
dt[, halfHour := hms::trunc_hms(hms, 30*60)]
# > plot LSOAs ----
p <- ggplot2::ggplot(dt,
aes(x = halfHour,
y = meankWh,
group = LSOA11CD,
colour = LSOA11CD)) +
geom_line() +
facet_grid(obsDate ~ .)
p <- p + theme(legend.position = "None") + # no legend
labs(x = "Time of Day",
y = "Mean household kWh per half hour",
caption = paste0("2011 Census & Observed Event Day kWh\n",
uniqueN(dt$LSOA11CD), " Southampton LSOAs")
)
ggplot2::ggsave(paste0(pPath,"lsoaProfiles_Event-Day-Sim.png"),
p)
return(dt)
}
# Function: SDRC 2.2 Baseline kWh winter weekday ----
# Need to calculate from weights file etc
getSDRC2.2BL <- function(){
sdrc2.2Path <- "~/Data/SAVE/processed/"
f <- "SAVEipfHouseholdWeights_Baseline_v_sdrc2.2_all_hhs.csv.gz" # get weights file
sdrc2.2_wts_OA_dt <- data.table::fread(paste0(sdrc2.2Path,
f))
# these are OA level weights so we need to aggregate (sum) them to LSOA level
sdrc2.2_wts_OA_dt[, zonecode := oaCode]
setkey(sdrc2.2_wts_OA_dt, zonecode)
summary(sdrc2.2_wts_OA_dt)
# add LSOA codes
sotonOA_11DT[, zonecode := OA11CD]
setkey(sotonOA_11DT, zonecode)
sdrc2.2_wts_OA_dt <- sdrc2.2_wts_OA_dt[sotonOA_11DT]
sdrc2.2_wts_LSOA_dt <- sdrc2.2_wts_OA_dt[ipfWeight > 0, .(wgt = sum(ipfWeight), # ignore zeros
mean_wgt = mean(ipfWeight),
sd_wgt = sd(ipfWeight)),
keyby = .(LSOA11CD,bmg_id)]
summary(sdrc2.2_wts_LSOA_dt)
# load the baseline January 2017 kWh data
f <- "jan2017clamp15minuteData.csv.gz" # get Wh data file
sdrc2.2_Jan2017kWh_dt <- data.table::fread(paste0(sdrc2.2Path,
f))
sdrc2.2_Jan2017kWh_dt[, kWh := consumptionWh/1000]
sdrc2.2_Jan2017kWh_dt[, obsDateTime := lubridate::as_datetime(obsDateTime,
tz = "UTC")] # make sure
# now create per hh means by half hour by weekday vs weekend for the whole period
# and then use this as basis for merge to weights
sdrc2.2_Jan2017kWh_dt[, wday := lubridate::wday(obsDateTime, label = TRUE)]
sdrc2.2_Jan2017kWh_dt[, weekday := "Weekday"]
sdrc2.2_Jan2017kWh_dt[wday == "Sat" | wday == "Sun", weekday := "Weekend"]
sdrc2.2_Jan2017kWh_dt[, hms := hms::as.hms(obsDateTime, tz = "UTC")] # crucial if you're doing the analysis outsde UK!
sdrc2.2_Jan2017kWh_dt[, halfHour := hms::trunc_hms(hms, 30*60)] # half hours
# sum to half-hours
sdrc2.2_Jan2017kWh_dt <- sdrc2.2_Jan2017kWh_dt[, .(nObs = .N,
sumkWh = sum(kWh, na.rm = TRUE)), # add up the two 15 min values
keyby = .(halfHour, obsDate = lubridate::as_date(obsDateTime),
weekday, bmg_id )]
# mean for weekday & weekend
aggJan2017kWh_dt <- sdrc2.2_Jan2017kWh_dt[, .(nObs = .N,
meankWh = mean(sumkWh, na.rm = TRUE)), # add up the two 15 min values
keyby = .(halfHour, weekday, bmg_id )]
# test plot of weekend vs weekday profiles
summary(aggJan2017kWh_dt)
# remove the NAs as they are when a bmg_id starts so we have no consumption history
plotDT <- sdrc2.2_Jan2017kWh_dt[, .(meankWh = mean(sumkWh)),
keyby = .(halfHour, weekday)]
p <- ggplot2::ggplot(plotDT, aes(x = halfHour, y = meankWh,
colour = weekday, group = weekday)) +
geom_line()
ggplot2::ggsave(paste0(pPath, "saveProfiles_Winter-Obs-kWh.png"), p)
setkey(aggJan2017kWh_dt, bmg_id)
setkey(sdrc2.2_wts_LSOA_dt, bmg_id)
aggJan2017kWh_dt[,s_halfhour := as.character(halfHour)] # so srvy doesn't break
sdrc2.2_wts_LSOA_dt <- sdrc2.2_wts_LSOA_dt[wgt > 0] # remove 0s
# merge kWh and weights
janSimdt <- aggJan2017kWh_dt[sdrc2.2_wts_LSOA_dt, allow.cartesian=TRUE]
# keep only what we need for weekdays
janSimWDaydt <- janSimdt[weekday == "Weekday", .(bmg_id, wgt, LSOA11CD,s_halfhour, weekday, meankWh)]
janSimWDaydt$weekday <- NULL # save some memory
# make a svydes
svy.wd_dt <- survey::svydesign(ids = ~bmg_id,
weights = ~wgt,
data = janSimWDaydt)
svyt_wd <- survey::svyby(~meankWh, ~s_halfhour+LSOA11CD, svy.wd_dt, # this is where you'd put the survey categories to average over
svymean,
keep.var = TRUE,
exclude=NULL,
na.action=na.pass
)
dt <- as.data.table(svyt_wd) # convert to data.table
dt[, halfHour := hms::parse_hms(s_halfhour)]
dt[, halfHour := hms::as.hms(halfHour, tzone = "UTC")]
# > plot LSOAs ----
p <- ggplot2::ggplot(dt, aes(x = halfHour,
y = meankWh, colour = LSOA11CD)) +
geom_line() +
theme(legend.position = "None")
ggplot2::ggsave(paste0(pPath,"lsoaProfiles_Winter-Weekdays-Obs-kWh-Sim.png"),
p)
return(dt)
}
# drake plan ----
plan <- drake::drake_plan(
sdrc2.1v1 = getSDRC2.1("v1.0"),
sdrc2.1v2 = getSDRC2.1("v2.0"),
sdrc2.2_ED = getSDRC2.2ED(), # loads from OA results file
#sdrc2.2_BL = getSDRC2.2BL(), # requires 2 inputs - weights & kWh obs
report = rmarkdown::render(
knitr_in(paste0(rmdRoot, "compareModels_IMA2019.Rmd")),
output_file = file_out(paste0(rmdRoot, "compareModels_IMA2019.html")),
quiet = FALSE
)
)
# test it ----
plan
config <- drake_config(plan)
vis_drake_graph(config)
# do it ----
make(plan)
---
title: "Comparing SAVE and GREEN Grid small area electricity demand models"
subtitle: "Analysis for paper presented at the IMA 7th World Congress, Galway June 19-21"
author: "Ben Anderson (University of Otago)"
date: 'Last run at: `r Sys.time()`'
output:
bookdown::html_document2:
code_folding: hide
fig_caption: yes
number_sections: yes
self_contained: no
toc: yes
toc_depth: 3
toc_float: yes
bibliography: '`r path.expand("~/bibliography.bib")`'
---
```{r setup, include=FALSE}
# This code should be called from the .R file in this directory to make sure all drake objects have already been created.
# change default
knitr::opts_chunk$set(echo = TRUE)
# Log compile time:
startTime <- proc.time()
# load packages we didn't load in the .R file
myPackages <- c(
"ggplot2", # plots
"kableExtra", # pretty tables
"skimr", # for skim
"viridis" # colour-blind friendly palettes
)
spatialec::loadLibraries(myPackages) # this will try to install packages it can't load
# ---- Local Settings ----
```
# Report Purpose
To compare SAVE small area models:
* SDRC2.1 v1.0: Time Use Model with all hot water = elec but not heat
* SDRC2.1 v2.0: Time Use Model with no hot water = elec and not heat
* SDRC2.2: Observed kWh Model
# Join models
Join SAVE models
```{r joinModels}
# Compare models ----
# join them at LSOA level
# v2.1 = 2001 v2,2 = 2011: some LSOAs may have changed
sdrc2.1v1_LSOA_dt <- drake::readd(sdrc2.1v1)
sdrc2.1v2_LSOA_dt <- drake::readd(sdrc2.1v2)
sdrc2.2_LSOA_dt <- drake::readd(sdrc2.2_BL)
sdrc2.1v1_LSOA_dt[, halfHour := hms] # mean weekday
sdrc2.1v1_LSOA_dt[, sdrc2.1v1meankWh := mean_powerkWh]
sdrc2.1v1_LSOA_dt <- sdrc2.1v1_LSOA_dt[,.(zonecode, halfHour, sdrc2.1v1meankWh)]
sdrc2.1v1_LSOA_dt[, halfHour := hms::trunc_hms(halfHour, 60*30)]
data.table::setkey(sdrc2.1v1_LSOA_dt, zonecode, halfHour)
sdrc2.1v2_LSOA_dt[, halfHour := hms] # mean weekday
sdrc2.1v2_LSOA_dt[, sdrc2.1v2meankWh := mean_powerkWh]
sdrc2.1v2_LSOA_dt <- sdrc2.1v2_LSOA_dt[,.(zonecode, halfHour, sdrc2.1v2meankWh)]
sdrc2.1v2_LSOA_dt[, halfHour := hms::trunc_hms(halfHour, 60*30)] # to remove the 1 sec (why is that there?)
data.table::setkey(sdrc2.1v2_LSOA_dt, zonecode, halfHour)
sdrc2.2_LSOA_dt[, sdrc2.2meankWh := meankWh] # mean weekday
sdrc2.2_LSOA_dt[, zonecode := LSOA11CD]
sdrc2.2_LSOA_dt <- sdrc2.2_LSOA_dt[,.(zonecode, halfHour, sdrc2.2meankWh)]
data.table::setkey(sdrc2.2_LSOA_dt, zonecode, halfHour)
# 2001 LSOAs:
data.table::uniqueN(sdrc2.1v1_LSOA_dt$zonecode)
data.table::uniqueN(sdrc2.1v2_LSOA_dt$zonecode)
# 2011 LSOAs:
data.table::uniqueN(sdrc2.2_LSOA_dt$zonecode)
# so we'll get some NAs through not matching
mergedLSOA <- sdrc2.1v1_LSOA_dt[sdrc2.1v2_LSOA_dt][sdrc2.2_LSOA_dt]
kableExtra::kable(head(mergedLSOA), caption = "Combined models data table: 10 rows") %>%
kable_styling()
t <- summary(mergedLSOA)
kableExtra::kable(t, caption = "Summary of combined models data table") %>%
kable_styling()
```
# Compare SDRC2.1 Model 1 with Model 2
At LSOA level...
```{r compare_v1_v2, fig.cap="Half-hourly mean kWh for model 1 vs model 2 for each LSOA"}
mergedLSOA[,peakFlag := "4. All other times"]
mergedLSOA[halfHour >= hms::as.hms("07:00:00") & halfHour < hms::as.hms("09:00:00"),
peakFlag := "1. Morning Peak"]
mergedLSOA[halfHour >= hms::as.hms("16:00:00") & halfHour < hms::as.hms("20:00:00"),
peakFlag := "3. Evening Peak"]
mergedLSOA[halfHour >= hms::as.hms("09:00:00") & halfHour <= hms::as.hms("16:00:00"),
peakFlag := "2. Day time"]
# plot LSOAs ----
p <- ggplot2::ggplot(mergedLSOA, aes(x = sdrc2.1v1meankWh,
y = sdrc2.1v2meankWh,
colour = zonecode)) +
geom_abline(intercept = 0, slope = 1) +
geom_point() +
facet_wrap(. ~ peakFlag) +
theme(legend.position = "None") +
labs(x = "mean kWh (Sim: Time Use model v1)",
y = "mean kWh (Sim: Time Use model v2)",
caption = paste0("Southampton: LSOAs\n",
"1:1 line shown"))
#xlim(0,0.8) +
#ylim(0,0.6) +
p
ggplot2::ggsave(paste0(pPath,"lsoaProfiles_Winter-Weekdays-Compare-SDCR2.1_v1_v2.png"), p)
peaks <- unique(mergedLSOA$peakFlag)
for(p in peaks){
message("Testing correlation for: ", p)
print(with(mergedLSOA[peakFlag == p], cor.test(sdrc2.1v1meankWh,
sdrc2.1v2meankWh,
na.action = "na.omit")
))
print(with(mergedLSOA[peakFlag == p], cor.test(sdrc2.1v1meankWh,
sdrc2.1v2meankWh,
na.action = "na.omit",
method = "spearman")
))
}
```
# Compare SDRC2.1 Model 1 with SDRC2.2 model
At LSOA level...
```{r compare_v1, fig.cap="Half-hourly mean kWh for model 1 vs observed for each LSOA"}
# plot LSOAs ----
p <- ggplot2::ggplot(mergedLSOA, aes(x = sdrc2.1v1meankWh,
y = sdrc2.2meankWh,
colour = zonecode)) +
geom_abline(intercept = 0, slope = 1) +
geom_point() +
facet_wrap(. ~ peakFlag) +
theme(legend.position = "None") +
labs(x = "mean kWh (Sim: Time Use model v1)",
y = "mean kWh (Sim: Observed kWh model)",
caption = paste0("Southampton: LSOAs\n",
"1:1 line shown"))
#xlim(0,0.8) +
#ylim(0,0.6) +
p
ggplot2::ggsave(paste0(pPath,"lsoaProfiles_Winter-Weekdays-Compare-SDCR2.1_v1_SDRC2.2.png"), p)
peaks <- unique(mergedLSOA$peakFlag)
for(p in peaks){
message("Testing correlation for: ", p)
print(with(mergedLSOA[peakFlag == p], cor.test(sdrc2.1v1meankWh,
sdrc2.2meankWh,
na.action = "na.omit")
))
print(with(mergedLSOA[peakFlag == p], cor.test(sdrc2.1v1meankWh,
sdrc2.2meankWh,
na.action = "na.omit",
method = "spearman")
))
}
```
# Compare SDRC2.1 Model 2 with SDRC2.2 model
At LSOA level...
```{r compare_v2, fig.cap="Half-hourly mean kWh for model 2 vs observed for each LSOA"}
# plot LSOAs ----
p <- ggplot2::ggplot(mergedLSOA, aes(x = sdrc2.1v2meankWh,
y = sdrc2.2meankWh,
colour = zonecode)) +
geom_abline(intercept = 0, slope = 1) +
geom_point() +
facet_wrap(. ~ peakFlag) +
theme(legend.position = "None") +
labs(x = "mean kWh (Sim: Time Use model v2)",
y = "mean kWh (Sim: Observed kWh model)",
caption = paste0("Southampton: LSOAs\n",
"1:1 line shown"))
#xlim(0,0.8) +
#ylim(0,0.6) +
p
ggplot2::ggsave(paste0(pPath,"lsoaProfiles_Winter-Weekdays-Compare-SDCR2.1_v2_SDRC2.2.png"), p)
peaks <- unique(mergedLSOA$peakFlag)
for(p in peaks){
message("Testing correlation for: ", p)
print(with(mergedLSOA[peakFlag == p], cor.test(sdrc2.1v2meankWh,
sdrc2.2meankWh,
na.action = "na.omit")
))
print(with(mergedLSOA[peakFlag == p], cor.test(sdrc2.1v2meankWh,
sdrc2.2meankWh,
na.action = "na.omit",
method = "spearman")
))
}
```
# Compare differences
```{r compare_delta, fig.cap="Half-hourly mean kWh for model 2 vs observed for each LSOA"}
mergedLSOA[, diffv1 := sdrc2.1v1meankWh - sdrc2.2meankWh]
mergedLSOA[, diffv2 := sdrc2.1v2meankWh - sdrc2.2meankWh]
t <- mergedLSOA[, .(meanDelta_v1 = mean(diffv1, na.rm = TRUE),
meanDelta_v2 = mean(diffv2, na.rm = TRUE)), keyby = .(peakFlag)]
kableExtra::kable(t, caption = "Summary of delta (+ve value means TU model higher)") %>%
kable_styling()
t1 <- mergedLSOA[, .(delta = diffv1, halfHour)]
t1[, label := "TU model 1: observed model"]
t2 <- mergedLSOA[, .(delta = diffv2, halfHour)]
t2[, label := "TU model 2 : observed model"]
dt <- rbind(t1,t2)
# plot LSOAs ----
p <- ggplot2::ggplot(dt, aes(x = halfHour, group = halfHour)) +
geom_boxplot(aes(y = delta)) +
facet_grid(label ~ .) +
labs(y = "delta mean kWh",
caption = paste0("Southampton: LSOAs"))
p
ggplot2::ggsave(paste0(pPath,"lsoaProfiles_Winter-Weekdays-Compare-delta_SDCR2.1_SDRC2.2.png"), p)
```
# Input data description
```{r v2.1v1}
skimr::skim(sdrc2.1v1_LSOA_dt)
```
```{r v2.1v2}
skimr::skim(sdrc2.1v2_LSOA_dt)
```
```{r v2.2}
skimr::skim(sdrc2.2_LSOA_dt)
```
# References
\ No newline at end of file
This diff is collapsed.
File deleted
File deleted
File deleted
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment