In CfSOtago/GREENGridEECA: Supports Analysis of NZ GREEN Grid data for EECA

knitr::opts_chunk$set(echo = FALSE) # by default turn off code echo

# Set start time ----
startTime <- proc.time()

# Local parameters ----
b2Kb <- 1024 #http://whatsabyte.com/P1/byteconverter.htm
b2Mb <- 1048576
plotLoc <- paste0(repoParams$repoLoc, "/docs/plots/") # where to put the plots (if any)

# Packages used in the report ----
rmdLibs <- c("ggplot2", # plots
          "kableExtra", # fancy tables
          "hms", # times
          "skimr", # for skim
          "tidyr", # for gather
          "dplyr", # for filter
          "readr"#for reading .csv
          )
GREENGridEECA::loadLibraries(rmdLibs)

# Local functions ----

\newpage

About

Citation

Report circulation:

Public – this report is intended for publication following EECA approval.

License

History

this report's edit history

Support

\newpage

Introduction

This report uses the GREEN Grid project [@stephenson_smart_2017] research data to analyse a variety of residential household appliances and their contribution to peak demand under several scenarios.

Data

# put all this here before you report on it
# this code _could_ go in the makeFile if you want...

#New DT
powerDT <- copy(origPowerDT) # need to copy as data.table works by reference https://www.rdocumentation.org/packages/data.table/versions/1.12.2/topics/copy

#DT modifications

#Excluding households based on Part A report: https://cfsotago.github.io/GREENGridEECA/#part-a--data-processing
exclude <- c("rf_14", "rf_25", "rf_26", "rf_43", "rf_46")
powerDT <- powerDT[!(linkID %in% exclude)]

# setting negative values to NA
powerDT <- powerDT[, meanPowerW := ifelse(meanPowerW <0, NA, meanPowerW)]
powerDT <- powerDT[, sdPowerW := ifelse(sdPowerW <0, NA, sdPowerW)]
powerDT <- powerDT[, minPowerW := ifelse(minPowerW <0, NA, minPowerW)]
powerDT <- powerDT[, maxPowerW := ifelse(maxPowerW <0, NA, maxPowerW)]


# set to NZ time from UTC
powerDT <- powerDT[, r_dateTime_nz := lubridate::as_datetime(r_dateTimeHalfHour, 
                                               tz = "Pacific/Auckland")] # this will be UTC unless you set this

#Define Winter/else
powerDT <- powerDT[, date := lubridate::date(r_dateTime_nz)]
powerDT <- powerDT[, obsHalfHour := hms::as.hms(r_dateTime_nz)]
#powerDT[, obsHalfHour := format(ymd_hms(r_dateTimeHalfHour), "%H:%M:%S")]
powerDT <- powerDT[, month := lubridate::month(r_dateTime_nz)]
powerDT <- powerDT[, peak := 0]

powerDT <- powerDT[month == 12 | month == 1 | month == 2, season := "Summer"]
powerDT <- powerDT[month == 3 | month == 4 | month == 5, season := "Autumn"]
powerDT <- powerDT[month == 6 | month == 7 | month == 8, season := "Winter"]
powerDT <- powerDT[month == 9 | month == 10 | month == 11, season := "Spring"]

#Setting times of peak demand 
OP1S <- hms::as.hms("00:00:00")
OP1E <- hms::as.hms("16:30:00")

PS <- hms::as.hms("17:00:00")
PE <- hms::as.hms("21:00:00")

OP2S <- hms::as.hms("21:30:00")
OP2E <- hms::as.hms("23:30:00")

powerDT <- powerDT[, peak := ifelse(obsHalfHour >= OP1S & obsHalfHour <= OP1E,
                                    "Off Peak 1",
                                    NA)]
powerDT <- powerDT[, peak := ifelse(obsHalfHour >= OP2S & obsHalfHour <= OP2E,
                                    "Off Peak 2",
                                    peak)]
powerDT <- powerDT[, peak := ifelse(obsHalfHour >= PS & obsHalfHour <= PE,
                                    "Peak",
                                    peak)]
# this stops the RHS coercion errors

Peak Contribution

For each dwelling using data averaged over 30 minute intervals to:

calculate total mean demand (kW) over the period
find the maximum total demand load (kW) during this period
calculate mean percentage contribution of each separately identifiable load type (lighting, hot water, heat pump, oven) in the maximum total load half-hour

We consider three different definition of peak:

Winter evenings 17:00 - 21:00
Regional/network peak
Sample co-incident peak - - select the 'peak' half-hour across all dwellings

Then for each dwelling: 1. calculate the mean demand over the peak period (if relevant) 2. find the maximum demand (kW) for the defined peak period 3. calculate percentage contribution

To do this we classify the circuit labels according to the load types mentioned above (see Table \@ref(tab:categoriseCircuitsTable) in the Data Annex).

Winter (June - August) evenings 17:00 – 21:00

In this section we select the peak-period half hours for all households for the winter evening period and then locate the half-hour which had the maximum total load for that dwelling. Table \@ref(tab:findMaxHalfHours) shows an example data extract.

powerDT[, year := lubridate::year(r_dateTime_nz)] # so we can do per-year analysis

selectDT <- powerDT[circuit == 'imputedTotalDemand_circuitsToSum_v1.1' & 
                   obsHalfHour >= hms::as.hms("17:00:00") & 
                     obsHalfHour <= hms::as.hms("21:00:00"), .SD[which.max(meanPowerW)], # finds only the rows which are the max of the mean total load (XX should we be using the maxPowerW?? XX)
                   by = .(linkID, season, year)] # max value within household, season, year

# This is the time of maximum peak demand for each household

#kableExtra::kable(head(selectDT), caption = "Extracted peak load half-hours - first 6 rows")  %>% # test result
#  kable_styling()

dataBucketDT <- data.table::data.table()
for(n in 1:nrow(selectDT)){
  # extract the rows which match the dateTime
  row <- selectDT[n]
  dt <- powerDT[linkID == row$linkID & # matches household
                  r_dateTime_nz == row$r_dateTime_nz, # matches dateTime
                .(linkID, year, season, r_dateTime_nz, circuit, nObs, meanPowerW, minPowerW, maxPowerW)] # selects just the basic columns we need
  dataBucketDT <- rbind(dataBucketDT, dt)
}



# extract _just_ the circuits which are (or contain):
# - total load
# - heat pumps
# - hot water
# - oven
# - lighting


dataBucketDT <- GREENGridEECA::labelEECACircuits(dataBucketDT) # use the function to split the circuit string for easier string matching

<<<<<<< HEAD
kableExtra::kable(head(dataBucketDT, 10), caption = "Example date peak load half-hour circuits - first 10 rows showing circuit demand at the peak total demand half-hour")  %>% # test result
=======
# test
t <- table(dataBucketDT$circuitLabel, dataBucketDT$eecaCircuit)

kableExtra::kable(t, caption = "Categorised circuits - XX check carefully! XX")  %>% # test result

kableExtra::kable(head(dataBucketDT, 10), caption = "Extracted peak load half-hour circuits - first 10 rows showing circuit demand at the peak total demand half-hour")  %>% # test result
>>>>>>> 3e35a9df1236bc2b754b9871020996aca93c5b03
  kable_styling()

# Full EECA category label check table is in Data Annex
categoriseCircuitsTable <- table(dataBucketDT$circuitLabel, dataBucketDT$eecaCircuit)

Table \@ref(tab:winterPeaktestTable) shows the winter peak load half-hours for each dwelling and year in winter. In this data extract, NaN indicates that no circuit data was obtained for this dwelling in this specific time period. We furthermore identify the contribution of each circuit label to the total demand. The complete data table is saved as "Extracted winter peak load half-hour circuits.csv" and is supplied with the report.

resAllSeasonsDT <- dcast(dataBucketDT,
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW") # take 1 household and turn data round

# Now just calculate the % using the identified circuits?


resAllSeasonsDT <- resAllSeasonsDT[, contribution_HW := `Hot water`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_HP := `Heat Pump or Heating`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_LI := `Lighting`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_OV := `Oven`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_XX := 1-((`Hot water` + `Heat Pump or Heating` + `Lighting` + `Oven`)/Total)]

kableExtra::kable(head(resAllSeasonsDT[season == "Winter"], 10), digits = 3,
                  caption = "Example data: Extracted winter peak load half-hour circuits (NaN indicates no such circuit)")  %>% # test result
  kable_styling()

write_csv(resAllSeasonsDT, "/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/csv_full_data_tables/Extracted winter peak load half-hour circuits.csv")

summaryP1DT <- copy(resAllSeasonsDT) #This is the summary table to present to EECA for the first definition of peak demand

summaryP1DT<-summaryP1DT[, c("Heat Pump or Heating", "Hot water", "Lighting", "Oven", "Total", "XX_Other") := NULL]

#summaryP1DT<-summaryP1DT[, averageContriHW := ]
summaryP1DT <- setnames(summaryP1DT, old=c("contribution_HW","contribution_HP", "contribution_LI", "contribution_OV", "contribution_XX"),
                        new=c("Hot Water", "Heat Pump", "Lighting", "Oven", "Others"))

#summaryP1DT <- data.table::as.data.table(gather(summaryP1DT, eecaCircuit, meanPowerW, `Heat Pump or Heating`:`XX_Other`, factor_key = TRUE))
summaryP1DT <- data.table::as.data.table(gather(summaryP1DT, eecaLabelContri, contribution, `Hot Water`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

summaryP1DT <- summaryP1DT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year)]

summaryP1DT <- summaryP1DT[, averageContri := averageContri*100 ]#Conversion to % for plot

summaryP1DT$season <- factor(summaryP1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot

#XX_eval =FALSE_XX

myPlot <- ggplot2::ggplot(summaryP1DT[season == "Winter" & year == 2015], aes(x= eecaLabelContri, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"), legend.position = "none",
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black")) +

  #facet_grid(season ~ .) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25, 30, 35)) +
  ggtitle("Contribution to peak demand (17:00-21:00) for winter in 2015") +
  labs(x='Circuit', y='Contribution in %') 



#XX_Into Annex??_XX

#XX_eval=FALSE_XX

#The stacks do not add up to 100 as we have averaged the contributions over all households. XXMention in the textXX

myPlot <- ggplot2::ggplot(summaryP1DT[year == 2015], aes(x= season, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_blank())+

    geom_text(aes(y = averageContri,
              label = sprintf("%0.2f",averageContri)),
           position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Season', y='Contribution in %') 
myPlot
#Now do the same for wattage contribution...

#This section aims to calculate the wattage contribution. We recycle previous results...

ContriW1DT <- dcast(dataBucketDT,#This tis the contribution DT for wattage and the first definition of peak demand
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW")
ContriW1DT<-ContriW1DT[, c( "Total") := NULL]#We do not need this column

ContriW1DT <- setnames(ContriW1DT, old=c("Heat Pump or Heating","Hot water", "Lighting", "Oven", "XX_Other"),
                        new=c("Heat Pump", "Hot Water", "Lighting", "Oven", "Others"))
ContriW1DT <- data.table::as.data.table(gather(ContriW1DT, eecaLabelContri, contribution, `Heat Pump`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

ContriW1DT <- ContriW1DT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year)]

ContriW1DT$season <- factor(ContriW1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriW1DT$eecaLabelContri <- factor(ContriW1DT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot

#Adding percentages

ContriW1DT <- ContriW1DT[year==2015]

ContriW1DT <- ContriW1DT[, totalW := sum(averageContri), keyby =.(season)]
ContriW1DT <- ContriW1DT[, percentage := (averageContri/totalW)*100]

ContriW1DT$season <- factor(ContriW1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriW1DT$eecaLabelContri <- factor(ContriW1DT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot

Figure \@ref(fig:wattageContributionPeak1Plot) and \@ref(fig:wattageContributionPeak1PlotPercent) show the contribution to peak demand of various appliances averaged over households. The contribution of hot water slightly increases in summer whereas contribution of heat pumps decreases over summer. Hot water is identified as the appliance that is relatively contributing the most to peak demand in 2015.

myPlot <- ggplot2::ggplot(ContriW1DT[year == 2015], aes(x= season, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_blank())+

    geom_text(aes(y = averageContri,
              label = sprintf("%0.2f",averageContri)),
           position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Season', y='Power in W/house') 
myPlot

myPlot <- ggplot2::ggplot(ContriW1DT[year == 2015], aes(x= "", y= percentage, fill = eecaLabelContri)) +
  geom_bar(stat = "identity", width = 1, color ="white")+
  #labs(x='', y='percentage')+
   facet_grid(facets=. ~season) +
  coord_polar("y", start = 0)+
  theme_void()+
 theme(text = element_text(family = "Cambria"),
       legend.position = 'bottom')+
  guides(fill=guide_legend(title="Circuit"))+
  geom_text(aes(y = percentage,
              label = sprintf("%0.0f%%",percentage)),
           position = position_stack(vjust = .5),
           size = 3)

myPlot

Region/network coincident peak

In this section we focus on two particular regions to analyse the contribution to peak deamand. These two regions are Hawke's Bay and Taranaki. Most of the household in our data set were monitored in these regions. In a first step, we link the region to the dwelling ID. We extract the power data for both regions seperately and exclude households that do not match with the defined regions.

regpeakDT <- copy(powerDT)

regpeakDT <- GREENGridEECA::labelEECACircuits(regpeakDT) # use the function to split the circuit string for easier string matching

regpeakDT <- regpeakDT[, circuitLabel := NULL]
regpeakDT <- regpeakDT[, peak := NULL]
regpeakDT <- regpeakDT[, month := NULL]
regpeakDT <- regpeakDT[, r_dateTimeHalfHour := NULL]
regpeakDT <- regpeakDT[, sdPowerW := NULL]
regpeakDT <- regpeakDT[, minPowerW := NULL]
regpeakDT <- regpeakDT[, circuitID := NULL]
regpeakDT <- regpeakDT[, circuit := NULL]
regpeakDT <- regpeakDT[, maxPowerW := NULL]
regpeakDT <- regpeakDT[, date := NULL]
regpeakDT <- regpeakDT[, obsHalfHour := NULL]#We do not need these, it slows down our computer...


regpeakDT <- regpeakDT[, region := NA]#Creating new variable

regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_08", "Taranaki", region)]#Match regions with IDs
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_26", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_21", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_13", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_24", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_27", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_19", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_11", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_10", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_22", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_06", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_14", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_23", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_07", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_25", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_12", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_16", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_20", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_09", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_15", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_17", "Taranaki", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_18", "Taranaki", region)]

regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_43", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_34", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_30", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_28", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_38", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_37", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_41", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_42", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_36", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_45", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_40", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_32", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_33", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_39", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_44", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_35", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_46", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_29", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_31", "Hawkes Bay", region)]
regpeakDT <- regpeakDT[, region := ifelse(`linkID`== "rf_47", "Hawkes Bay", region)]


ind <- with(regpeakDT, region == is.na(regpeakDT$region))#Now let's remove the households that do not match the two regions
regpeakDT <- regpeakDT[!ind, ]#Removing households that do not match and have NA as value

taranakiDT <- copy(regpeakDT)
taranakiDT <- taranakiDT[region=="Taranaki"]#taranaki only
hawkesDT <- copy(regpeakDT)
hawkesDT <- hawkesDT[region=="Hawkes Bay"]#hawkes bay only

##Data import
GXPtaranakiDT <- data.table::as.data.table(
 readr::read_csv("/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/taranakiGxpTop100DT.csv",
                       col_types = cols(
                       region = col_character(),
                       regionName = col_character(),
                       rDateTime = col_datetime(format = ""),
                       weekdays = col_character(),
                       peakPeriod = col_character(),
                       sumkWh = col_double(),
                       nObs = col_double(),
                       nRegions = col_double(),
                       nPOCs = col_double(),
                       nNetworks = col_double(),
                       nGenTypes = col_double(),
                       month = col_character() # <- crucial otherwise readr infers integers for some reason
                     )
                 )
)

# set to NZ time from UTC
GXPtaranakiDT <- GXPtaranakiDT[, r_dateTime_nz := lubridate::as_datetime(rDateTime, 
                                               tz = "Pacific/Auckland")] # this will be UTC unless you set this

#Define Winter/else
GXPtaranakiDT <- GXPtaranakiDT[, date := lubridate::date(r_dateTime_nz)]
GXPtaranakiDT <- GXPtaranakiDT[, obsHalfHour := hms::as.hms(r_dateTime_nz)]
#powerDT[, obsHalfHour := format(ymd_hms(r_dateTimeHalfHour), "%H:%M:%S")]
GXPtaranakiDT <- GXPtaranakiDT[, month := lubridate::month(r_dateTime_nz)]


GXPtaranakiDT <- GXPtaranakiDT[month == 12 | month == 1 | month == 2, season := "Summer"]
GXPtaranakiDT <- GXPtaranakiDT[month == 3 | month == 4 | month == 5, season := "Autumn"]
GXPtaranakiDT <- GXPtaranakiDT[month == 6 | month == 7 | month == 8, season := "Winter"]
GXPtaranakiDT <- GXPtaranakiDT[month == 9 | month == 10 | month == 11, season := "Spring"]

##Data import
GXPHawkesDT <- data.table::as.data.table(
 readr::read_csv("/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/hawkesBayGxpTop100DT.csv",
                       col_types = cols(
                       region = col_character(),
                       regionName = col_character(),
                       rDateTime = col_datetime(format = ""),
                       weekdays = col_character(),
                       peakPeriod = col_character(),
                       sumkWh = col_double(),
                       nObs = col_double(),
                       nRegions = col_double(),
                       nPOCs = col_double(),
                       nNetworks = col_double(),
                       nGenTypes = col_double(),
                       month = col_character() # <- crucial otherwise readr infers integers for some reason
                     )
                 )
)

# set to NZ time from UTC
GXPHawkesDT <- GXPHawkesDT[, r_dateTime_nz := lubridate::as_datetime(rDateTime, 
                                               tz = "Pacific/Auckland")] # this will be UTC unless you set this

#Define Winter/else
GXPHawkesDT <- GXPHawkesDT[, date := lubridate::date(r_dateTime_nz)]
GXPHawkesDT <- GXPHawkesDT[, obsHalfHour := hms::as.hms(r_dateTime_nz)]
#powerDT[, obsHalfHour := format(ymd_hms(r_dateTimeHalfHour), "%H:%M:%S")]
GXPHawkesDT <- GXPHawkesDT[, month := lubridate::month(r_dateTime_nz)]


GXPHawkesDT <- GXPHawkesDT[month == 12 | month == 1 | month == 2, season := "Summer"]
GXPHawkesDT <- GXPHawkesDT[month == 3 | month == 4 | month == 5, season := "Autumn"]
GXPHawkesDT <- GXPHawkesDT[month == 6 | month == 7 | month == 8, season := "Winter"]
GXPHawkesDT <- GXPHawkesDT[month == 9 | month == 10 | month == 11, season := "Spring"]

Data on regional demand for the Hawke's Bay and Taranaki Grid Exit Point (GXP) was obtained from the Electriicty Authority and the top 100 demand half hours (i.e. the date time of the top 100 half hours) were extracted from this data. The demand of each house in the two regions was then extracted for the top 100 half-hours.

#These loops take one night to process....
#I set eval=FALSE to prevent it from running
#They look at the GXP rDateTime and then find this time in the region DT and extract the whole row into a new DT

#Hawke's Bay
extractHawkesDT <- NULL
dateTimelistHawkes <- hawkesDT[["r_dateTime_nz"]]
dateTimelist <- GXPHawkesDT[["r_dateTime_nz"]]
i<-0
for (time in c(dateTimelist)){
  for (newtime in c(dateTimelistHawkes)){
    i<-i+1
    if (toString(time) == toString(newtime)){
      #test = do.call("rbind",(hawkesDT[i,]))
     extractHawkesDT = rbind(extractHawkesDT, data.table(hawkesDT[i,]))
     # print(c(hawkesDT[i,]))
    }
  }
  i <- 0
}

#Taranaki

extractTaranakiDT <- NULL
dateTimelistTaranaki <- taranakiDT[["r_dateTime_nz"]]
dateTimelist <- GXPtaranakiDT[["r_dateTime_nz"]]
i<-0
for (time in c(dateTimelist)){
  for (newtime in c(dateTimelistTaranaki)){
    i<-i+1
    if (toString(time) == toString(newtime)){
      #test = do.call("rbind",(hawkesDT[i,]))
      extractTaranakiDT = rbind(extractTaranakiDT, data.table(taranakiDT[i,]))
     # print(c(hawkesDT[i,]))
    }
  }
  i <- 0
}

#Maybe it is a good idea to export these two tables as a csv in case we want to re-run analysis without waiting for ages.

#write.csv(extractHawkesDT, file = "extractHawkes.csv", row.names = FALSE) #They are stored in the Dropbox data folder

#write.csv(extractTaranakiDT, file = "extractTaranaki.csv", row.names = FALSE)

#The loops above do not allow us to do doReport() in a reasonable time. I therefore decided to export the result and import it rather than calculating it again in the loops.

#Loading Hawkes Bay extract data
extractHawkesDT <- data.table::as.data.table(
 readr::read_csv("/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/Matched_date_times_top 100/extractHawkes.csv",
                       col_types = cols(
                       linkID = col_character(),
                       nObs = col_double(),
                       r_dateTime_nz = col_datetime(format = ""),
                       meanPowerW = col_double(),
                       season = col_character(),
                       year = col_character(),
                       eecaCircuit = col_character(),
                       region = col_character() # <- crucial otherwise readr infers integers for some reason
                     )
                 )
)


#Loading Taranaki extract data

extractTaranakiDT <- data.table::as.data.table(
 readr::read_csv("/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/Matched_date_times_top 100/extractTaranaki.csv",
                       col_types = cols(
                       linkID = col_character(),
                       nObs = col_double(),
                       r_dateTime_nz = col_datetime(format = ""),
                       meanPowerW = col_double(),
                       season = col_character(),
                       year = col_character(),
                       eecaCircuit = col_character(),
                       region = col_character() # <- crucial otherwise readr infers integers for some reason
                     )
                 )
)

Hawke's Bay

Table \@ref(tab:RegionalAnalysisHawkesContribution) presents example data of the successfully extracted peak load half-hours by circuit for 2015. Peak in the following is defined as the period of time where the to 100 GXP date times occur. The complete data table is saved as "Extracted peak load half-hour circuits for Hawke's Bay.csv" and is supplied with the report.

resAllSeasonsDT <- copy(extractHawkesDT)

resAllSeasonsDT <- dcast(resAllSeasonsDT,
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW") # take 1 household and turn data round

# Now just calculate the % using the identified circuits?


resAllSeasonsDT <- resAllSeasonsDT[, contribution_HW := `Hot water`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_HP := `Heat Pump or Heating`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_LI := `Lighting`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_OV := `Oven`/Total]
#resAllSeasonsDT <- resAllSeasonsDT[, contribution_XX := 1-((`Hot water` + `Heat Pump or Heating` + `Lighting` + `Oven`)/Total)]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_XX := `XX_Other`/Total]

kableExtra::kable(head(resAllSeasonsDT[year == "2015"], 10), digits = 3,
                  caption = "Example data: Extracted peak load half-hour circuits for Hawke's Bay (NaN indicates no such circuit)")  %>% # test result
  kable_styling()

write_csv(resAllSeasonsDT, "/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/csv_full_data_tables/Extracted peak load half-hour circuits for Hawke's Bay.csv")

#Eval=False, we do not need this 

summaryP1DT <- copy(resAllSeasonsDT) #This is the summary table to present to EECA for the first definition of peak demand

summaryP1DT<-summaryP1DT[, c("Heat Pump or Heating", "Hot water", "Lighting", "Oven", "Total", "XX_Other") := NULL]

#summaryP1DT<-summaryP1DT[, averageContriHW := ]
summaryP1DT <- setnames(summaryP1DT, old=c("contribution_HW","contribution_HP", "contribution_LI", "contribution_OV", "contribution_XX"),
                        new=c("Hot Water", "Heat Pump", "Lighting", "Oven", "Others"))

#summaryP1DT <- data.table::as.data.table(gather(summaryP1DT, eecaCircuit, meanPowerW, `Heat Pump or Heating`:`XX_Other`, factor_key = TRUE))
summaryP1DT <- data.table::as.data.table(gather(summaryP1DT, eecaLabelContri, contribution, `Hot Water`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

summaryP1DT <- summaryP1DT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year)]

summaryP1DT <- summaryP1DT[, averageContri := averageContri*100]#Conversion to % for plot

summaryP1DT$season <- factor(summaryP1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot

#Eval=false, we use the pie instead

#The stacks do not add up to 100 as we have averaged the contributions over all households. XXMention in the textXX XXDue to the fact that all top100 half hours were in winter we only have winter in the plotXX

myPlot <- ggplot2::ggplot(summaryP1DT[year == 2015], aes(x= season, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_blank())+

    geom_text(aes(y = averageContri,
              label = sprintf("%0.2f",averageContri)),
           position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Season', y='Contribution in %') 
myPlot
#Now do the same for wattage contribution...

#This section aims to calculate the wattage contribution. We recycle previous results...
ContriW1DT <- copy(extractHawkesDT)

ContriW1DT <- dcast(ContriW1DT,#This tis the contribution DT for wattage and the first definition of peak demand
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW")
ContriW1DT<-ContriW1DT[, c( "Total") := NULL]#We do not need this column

ContriW1DT <- setnames(ContriW1DT, old=c("Heat Pump or Heating","Hot water", "Lighting", "Oven", "XX_Other"),
                        new=c("Heat Pump", "Hot Water", "Lighting", "Oven", "Others"))
ContriW1DT <- data.table::as.data.table(gather(ContriW1DT, eecaLabelContri, contribution, `Heat Pump`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

ContriW1DT <- ContriW1DT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year)]

ContriW1DT$season <- factor(ContriW1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriW1DT$eecaLabelContri <- factor(ContriW1DT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot


#Adding percentages

ContriW1DT <- ContriW1DT[year==2015]

ContriW1DT <- ContriW1DT[, totalW := sum(averageContri), keyby =.(season)]
ContriW1DT <- ContriW1DT[, percentage := (averageContri/totalW)*100]

ContriW1DT$season <- factor(ContriW1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriW1DT$eecaLabelContri <- factor(ContriW1DT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot

In terms of average wattage contribution of the top 100 GXP date times in Hawke's Bay (Figure \@ref(fig:wattageContributionHawkesPlot) ), Hot Water and Heat Pumps contribute 682 Watts and 672 Watts respectively. Heat Pumps and Hot Water contributed almost equally to peak demand in Hawke's Bay (27% and 28%), followed by Lighitng (10%) and Oven (3%). Other demand that was not further specified contributed 32%. This is presented in Figure \@ref(fig:wattageContributionHawkesPlotPercentage).

myPlot <- ggplot2::ggplot(ContriW1DT[year == 2015], aes(x= season, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_blank())+

    geom_text(aes(y = averageContri,
              label = sprintf("%0.2f",averageContri)),
           position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Season', y='Power in W/house') 
myPlot

myPlot <- ggplot2::ggplot(ContriW1DT[year == 2015], aes(x= "", y= percentage, fill = eecaLabelContri)) +
  geom_bar(stat = "identity", width = 1, color ="white")+
  #labs(x='', y='percentage')+
   facet_grid(facets=. ~season) +
  coord_polar("y", start = 0)+
  theme_void()+
 theme(text = element_text(family = "Cambria"),
       legend.position = 'bottom')+
  guides(fill=guide_legend(title="Circuit"))+
  geom_text(aes(y = percentage,
              label = sprintf("%0.0f%%",percentage)),
           position = position_stack(vjust = .5),
           size = 4)

myPlot

Taranaki

Table \@ref(tab:RegionalAnalysisTaranakiContribution) presents example data of extracted peak load half-hours by circuit for 2015 in Taranaki. The complete data table is saved as "Extracted peak load half-hour circuits for Taranaki.csv" and is supplied with the report.

resAllSeasonsDT <- copy(extractTaranakiDT)

resAllSeasonsDT <- dcast(resAllSeasonsDT,
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW") # take 1 household and turn data round

# Now just calculate the % using the identified circuits?


resAllSeasonsDT <- resAllSeasonsDT[, contribution_HW := `Hot water`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_HP := `Heat Pump or Heating`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_LI := `Lighting`/Total]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_OV := `Oven`/Total]
#resAllSeasonsDT <- resAllSeasonsDT[, contribution_XX := 1-((`Hot water` + `Heat Pump or Heating` + `Lighting` + `Oven`)/Total)]
resAllSeasonsDT <- resAllSeasonsDT[, contribution_XX := `XX_Other`/Total]

kableExtra::kable(head(resAllSeasonsDT[year == "2015"], 10), digits = 3,
                  caption = "Example data: Extracted peak load half-hour circuits (NaN indicates no such circuit)")  %>% # test result
  kable_styling()


write_csv(resAllSeasonsDT, "/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/csv_full_data_tables/Extracted peak load half-hour circuits for Taranaki.csv")

#Eval=False, not needed
summaryP1DT <- copy(resAllSeasonsDT) #This is the summary table to present to EECA for the first definition of peak demand

summaryP1DT<-summaryP1DT[, c("Heat Pump or Heating", "Hot water", "Lighting", "Oven", "Total", "XX_Other") := NULL]

#summaryP1DT<-summaryP1DT[, averageContriHW := ]
summaryP1DT <- setnames(summaryP1DT, old=c("contribution_HW","contribution_HP", "contribution_LI", "contribution_OV", "contribution_XX"),
                        new=c("Hot Water", "Heat Pump", "Lighting", "Oven", "Others"))

#summaryP1DT <- data.table::as.data.table(gather(summaryP1DT, eecaCircuit, meanPowerW, `Heat Pump or Heating`:`XX_Other`, factor_key = TRUE))
summaryP1DT <- data.table::as.data.table(gather(summaryP1DT, eecaLabelContri, contribution, `Hot Water`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

summaryP1DT <- summaryP1DT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year)]

summaryP1DT <- summaryP1DT[, averageContri := averageContri*100]#Conversion to % for plot

summaryP1DT$season <- factor(summaryP1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot

#Adding percentages

summaryP1DT <- summaryP1DT[year==2015]

summaryP1DT <- summaryP1DT[, totalW := sum(averageContri), keyby =.(season)]
summaryP1DT <- summaryP1DT[, percentage := (averageContri/totalW)*100]

summaryP1DT$season <- factor(summaryP1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
summaryP1DT$eecaLabelContri <- factor(summaryP1DT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot

#Eval=False, we use the pie chart instead

#The stacks do not add up to 100 as we have averaged the contributions over all households. XXMention in the textXX XXDue to the fact that all top100 half hours were in winter we only have winter in the plotXX

myPlot <- ggplot2::ggplot(summaryP1DT[year == 2015], aes(x= season, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_blank())+

    geom_text(aes(y = averageContri,
              label = sprintf("%0.2f",averageContri)),
           position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Season', y='Contribution in %') 
myPlot
#Now do the same for wattage contribution...

#This section aims to calculate the wattage contribution. We recycle previous results...
ContriW1DT <- copy(extractTaranakiDT)

ContriW1DT <- dcast(ContriW1DT,#This tis the contribution DT for wattage and the first definition of peak demand
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW")
ContriW1DT<-ContriW1DT[, c( "Total") := NULL]#We do not need this column

ContriW1DT <- setnames(ContriW1DT, old=c("Heat Pump or Heating","Hot water", "Lighting", "Oven", "XX_Other"),
                        new=c("Heat Pump", "Hot Water", "Lighting", "Oven", "Others"))
ContriW1DT <- data.table::as.data.table(gather(ContriW1DT, eecaLabelContri, contribution, `Heat Pump`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

ContriW1DT <- ContriW1DT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year)]

ContriW1DT$season <- factor(ContriW1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriW1DT$eecaLabelContri <- factor(ContriW1DT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot



#Adding percentages

ContriW1DT <- ContriW1DT[year==2015]

ContriW1DT <- ContriW1DT[, totalW := sum(averageContri), keyby =.(season)]
ContriW1DT <- ContriW1DT[, percentage := (averageContri/totalW)*100]

ContriW1DT$season <- factor(ContriW1DT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriW1DT$eecaLabelContri <- factor(ContriW1DT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot

The averaged wattage contribution by circuit is presented in Figure \@ref(fig:wattageContributionTaranakiPlot).Figure \@ref(fig:wattageContributionTaranakiPlotPercentage) shows the percentage contribution to peak demand as an average over the date times in the GXP data by circuits. Lighting contributes the most to peak demand with 25%, followed by Heat Pumps (23%), Hot Water (14%), and Oven (8%). Not further specified loads contribute 30% to peak demand in Taranaki.

myPlot <- ggplot2::ggplot(ContriW1DT[year == 2015], aes(x= season, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_blank())+

    geom_text(aes(y = averageContri,
              label = sprintf("%0.2f",averageContri)),
           position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Season', y='Power in W/house') 
myPlot

myPlot <- ggplot2::ggplot(ContriW1DT[year == 2015], aes(x= "", y= percentage, fill = eecaLabelContri)) +
  geom_bar(stat = "identity", width = 1, color ="white")+
  #labs(x='', y='percentage')+
   facet_grid(facets=. ~season) +
  coord_polar("y", start = 0)+
  theme_void()+
 theme(text = element_text(family = "Cambria"),
       legend.position = 'bottom')+
  guides(fill=guide_legend(title="Circuit"))+
  geom_text(aes(y = percentage,
              label = sprintf("%0.0f%%",percentage)),
           position = position_stack(vjust = .5),
           size = 4)

myPlot

Sample coincident peak

The sample coincident peak utilizes the GREENGrid sample to define peak demand. This is in contrast to section 4.1 where peak demand was predefined as the time between 17:00 and 21:00 in winter and section 4.2 wehere peak demand was predefined as the top 100 half hours based on regional GXP data.

Table \@ref(tab:sampleCoincidentPeak) identifies the half hours of peak demand in 2015 for the total sample load. Results are presented in a descending order. Please note that the aim of this data table is to identify the date and time of overall sample peak demand. The complete data table is saved as "Extracted winter peak load half-hours in descending order for 2015" and is supplied with the report.

To do that, we have to consider the number of observations for each half hour (otherwise a high number of observations would have a higher impact on overall demand). The displayed average power in Watts is, therefore, not the total for the half-hour but a division of the sum of wattage per half-hour and the number of observations. We would like to calculate the relative demand that considers a varying number of observations (nobs).

sampleDT <- labelEECACircuits(powerDT)

sampleDT <- sampleDT[, summeanPowerW := NA]

sampleShortDT <- sampleDT[, .(summeanPowerW = sum(meanPowerW, na.rm = TRUE), 
                              nobs=.N), keyby=.(eecaCircuit, date, season, year, obsHalfHour)]#This is the aggregated Wattage for each circuit month etc. see keyby in order to conduct sample coincident peak analysis 


sampleShortDT <- sampleShortDT[, averagePowerW := summeanPowerW/nobs]#We have to consider that some of the obs have only one obs, some have up to 200. We need to take the average to identify the top 20 entries (their time of occurance)


sample2015DT <- sampleShortDT[year == 2015] # Getting rid of years that are not necessary for this analysis

#sample2015DT <- sort(sample2015DT$averagePowerW, decreasing = TRUE)

sample2015DT <- setorder(sample2015DT, -averagePowerW)

kableExtra::kable(head(sample2015DT[year == "2015"], 10), digits = 3,
                  caption = "Example data: Extracted winter peak load half-hours in descending order for 2015")  %>% # test result
  kable_styling()


write_csv(resAllSeasonsDT, "/Users/carsten.dortans/Dropbox/GreenGrid_EECA/data/csv_full_data_tables/Extracted winter peak load half-hours in descending order for 2015.csv")


#Now identufy these top20 times in the original data table and conduct analysis again. We only use this data table to identify the peak half hours. Now we go back to original data files

####Filtering the top20 peak half-hours based on our sample. XXThis is hard-coded - is there a better way to do it??XX

In summary, the following half hours were identified as the top 20 sample peak half hours in 2015:

date == "2015-06-28" & obsHalfHour == hms::as.hms("18:00:00")
date == "2015-08-21" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-06-23" & obsHalfHour == hms::as.hms("18:00:00")
date == "2015-06-24" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-05-29" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-06-17" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-08-12" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-06-23" & obsHalfHour == hms::as.hms("17:30:00")
date == "2015-06-26" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-06-25" & obsHalfHour == hms::as.hms("07:30:00")
date == "2015-06-23" & obsHalfHour == hms::as.hms("18:30:00")
date == "2015-05-25" & obsHalfHour == hms::as.hms("19:00:00")
date == "2015-06-23" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-08-10" & obsHalfHour == hms::as.hms("20:00:00")
date == "2015-06-26" & obsHalfHour == hms::as.hms("07:30:00")
date == "2015-07-10" & obsHalfHour == hms::as.hms("07:30:00")
date == "2015-08-28" & obsHalfHour == hms::as.hms("07:00:00")
date == "2015-07-09" & obsHalfHour == hms::as.hms("07:30:00")
date == "2015-06-23" & obsHalfHour == hms::as.hms("19:30:00")
date == "2015-05-27" & obsHalfHour == hms::as.hms("07:00:00")

filteredDT <- powerDT %>% filter(date == "2015-06-28" & obsHalfHour == hms::as.hms("18:00:00") |
                                   date == "2015-08-21" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-06-23" & obsHalfHour == hms::as.hms("18:00:00") |
                                   date == "2015-06-24" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-05-29" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-06-17" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-08-12" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-06-23" & obsHalfHour == hms::as.hms("17:30:00") |
                                   date == "2015-06-26" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-06-25" & obsHalfHour == hms::as.hms("07:30:00") |
                                   date == "2015-06-23" & obsHalfHour == hms::as.hms("18:30:00") |
                                   date == "2015-05-25" & obsHalfHour == hms::as.hms("19:00:00") |
                                   date == "2015-06-23" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-08-10" & obsHalfHour == hms::as.hms("20:00:00") |
                                   date == "2015-06-26" & obsHalfHour == hms::as.hms("07:30:00") |
                                   date == "2015-07-10" & obsHalfHour == hms::as.hms("07:30:00") |
                                   date == "2015-08-28" & obsHalfHour == hms::as.hms("07:00:00") |
                                   date == "2015-07-09" & obsHalfHour == hms::as.hms("07:30:00") |
                                   date == "2015-06-23" & obsHalfHour == hms::as.hms("19:30:00") |
                                   date == "2015-05-27" & obsHalfHour == hms::as.hms("07:00:00"))

filteredDT <- as.data.table(filteredDT)
filteredDT <- labelEECACircuits(filteredDT)

contriSampleDT <- dcast(filteredDT,
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW") # take 1 household and turn data round

# Now just calculate the % using the identified circuits?


contriSampleDT <- contriSampleDT[, contribution_HW := `Hot water`/Total]
contriSampleDT <- contriSampleDT[, contribution_HP := `Heat Pump or Heating`/Total]
contriSampleDT <- contriSampleDT[, contribution_LI := `Lighting`/Total]
contriSampleDT <- contriSampleDT[, contribution_OV := `Oven`/Total]
contriSampleDT <- contriSampleDT[, contribution_XX := 1-((`Hot water` + `Heat Pump or Heating` + `Lighting` + `Oven`)/Total)]

#Eval=False, we use the pie chart instead
summaryPSDT <- copy(contriSampleDT) #This is the summary table to present to EECA for the first definition of peak demand

summaryPSDT<-summaryPSDT[, c("Heat Pump or Heating", "Hot water", "Lighting", "Oven", "Total", "XX_Other") := NULL]


summaryPSDT <- setnames(summaryPSDT, old=c("contribution_HW","contribution_HP", "contribution_LI", "contribution_OV", "contribution_XX"),
                        new=c("Hot Water", "Heat Pump", "Lighting", "Oven", "Others"))

summaryPSDT <- data.table::as.data.table(gather(summaryPSDT, eecaLabelContri, contribution, `Hot Water`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

summaryPSDT <- summaryPSDT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year, r_dateTime_nz)]

summaryPSDT <- summaryPSDT[, averageContri := averageContri*100 ]#Conversion to % for plot

summaryPSDT$season <- factor(summaryPSDT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot

Figure \@ref(fig:wattageContributionPeakSamplePlot) presents the wattage contribution on the 2015-06-28 at 18:00:00 as an average over the households for one of the identified top peak half hours. The results show that Hot Water contributes the most, followed by Heat Pumps and Lighting. The plot also visualizes that much power was drawn from circuits other than the identified ones (summarized in "Others"). Figure \@ref(fig:wattageContributionPeakSamplePlotPercentage) shows the associated percentage contribution. Since the chosen half-hour of peak demand is on the 2015-06-28 at 18:00:00 (winter) Heat Pump and Lighting demand are clearly contributing to peak demand. Furthermore, How Water is the circuit contributing the most.

#This section aims to calculate the wattage contribution. We recycle previous results...

ContriWSDT <- dcast(filteredDT,#This tis the contribution DT for wattage and the first definition of peak demand
           linkID + year + season + r_dateTime_nz ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "meanPowerW")
ContriWSDT<-ContriWSDT[, c( "Total") := NULL]#We do not need this column

ContriWSDT <- setnames(ContriWSDT, old=c("Heat Pump or Heating","Hot water", "Lighting", "Oven", "XX_Other"),
                        new=c("Heat Pump", "Hot Water", "Lighting", "Oven", "Others"))
ContriWSDT <- data.table::as.data.table(gather(ContriWSDT, eecaLabelContri, contribution, `Heat Pump`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

ContriWSDT <- ContriWSDT[, .(averageContri = mean(contribution, na.rm = TRUE), 
                              nobs=.N, sd=sd(contribution, na.rm = TRUE)), keyby=.(eecaLabelContri, season, year, r_dateTime_nz)]

ContriWSDT$season <- factor(ContriWSDT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriWSDT$eecaLabelContri <- factor(ContriWSDT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot



#Adding percentages

ContriWSDT <- ContriWSDT[r_dateTime_nz == "2015-06-28 18:00:00"]

ContriWSDT <- ContriWSDT[, totalW := sum(averageContri), keyby =.(season)]
ContriWSDT <- ContriWSDT[, percentage := (averageContri/totalW)*100]

ContriWSDT$season <- factor(ContriWSDT$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot
ContriWSDT$eecaLabelContri <- factor(ContriWSDT$eecaLabelContri, levels = c("Hot Water","Heat Pump",
                                                    "Lighting", "Oven", "Others"))#This puts the seasons in the right order in the plot

myPlot <- ggplot2::ggplot(ContriWSDT[r_dateTime_nz == "2015-06-28 18:00:00"], aes(x= r_dateTime_nz, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_blank())+
    geom_text(aes(y = averageContri,
             label = sprintf("%0.2f",averageContri)),
            position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x=' Example date time: 2015-06-28 18:00:00 (one of the top 20 half hours)', y='Power in W/house') 
myPlot

myPlot <- ggplot2::ggplot(ContriWSDT[r_dateTime_nz == "2015-06-28 18:00:00"], aes(x= "", y= percentage, fill = eecaLabelContri)) +
  geom_bar(stat = "identity", width = 1, color ="white")+
  #labs(x='', y='percentage')+
   facet_grid(facets=. ~season) +
  coord_polar("y", start = 0)+
  theme_void()+
 theme(text = element_text(family = "Cambria"),
       legend.position = 'bottom')+
  guides(fill=guide_legend(title="Circuit"))+
  geom_text(aes(y = percentage,
              label = sprintf("%0.0f%%",percentage)),
           position = position_stack(vjust = .5),
           size = 4)

myPlot

#EVAL=FALSE We use the pie chart instead

#The stacks do not add up to 100 as we have averaged the contributions over all households. XXMention in the textXX

myPlot <- ggplot2::ggplot(summaryPSDT[r_dateTime_nz =="2015-06-28 18:00:00"], aes(x= eecaLabelContri, fill = eecaLabelContri)) +
  geom_bar(aes(y = averageContri), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        legend.position = "none",
        axis.ticks.y = element_blank())+

    geom_text(aes(y = averageContri,
             label = sprintf("%0.2f",averageContri)),
            position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Circuit', y='Contribution in %') 
myPlot
#Now do the same for wattage contribution...

Annual Contribution

This section presents annual and seasonal energy consumption for a variety of circuits and the variability between dwellings. This analysis considers the year 2015.

Figure \@ref(fig:annualplot) shows the distribution of annual energy consumption (kWh) between houses for various circuits. Hot Water draws the most energy per year with a pronounced variation between houses. Lighting, Heat Pump, and Oven follow. Not further specified circuits are summarised in "Others".

# kWh

# new DT ( do we need to?)
energyDT <- copy(powerDT)

#energyDT <- energyDT[, meanEnergyWh := meanPowerW/2]#Converting W into Wh energy 



#Conversion to energy Wh
energyDT <- energyDT[, meanEnergyWh := meanPowerW/2]
analysisDT <- copy(energyDT) # another new DT. Do we need to?
analysisDT <- analysisDT[, totalEnergyWh := sum(meanEnergyWh, na.rm = TRUE), 
                         keyby = .(linkID, eecaCircuit, year)]#Ignores NAs
analysisShortDT <- analysisDT[, .(totalEnergyWh = totalEnergyWh), 
                              keyby = .(linkID, eecaCircuit, year)] #For some reason by does not remove duplicates without crashing....

analysisShortDT <- unique(analysisShortDT)# Removing duplicates

# Calculating the contribution to annual energy consumption and turning off e-notation

options(scipen = 999)

analysisShortDT <- dcast(analysisShortDT, 
           linkID + year ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "totalEnergyWh") # take 1 household and turn data round

analysisShortDT <- analysisShortDT[, contribution_HW := `Hot water`/Total]
analysisShortDT <- analysisShortDT[, contribution_HP := `Heat Pump or Heating`/Total]
analysisShortDT <- analysisShortDT[, contribution_LI := `Lighting`/Total]
analysisShortDT <- analysisShortDT[, contribution_OV := `Oven`/Total]
analysisShortDT <- analysisShortDT[, contribution_XX := 1-((`Hot water` + `Heat Pump or Heating` + `Lighting` + `Oven`)/Total)]
# There are a lot of NA/Nan in the data. XXSensible to turn NaN into NA to receive more results??XX

analysisShortDT <- setnames(analysisShortDT, old=c("Heat Pump or Heating","Hot water", "Lighting", "Oven", "XX_Other"),
                        new=c("Heat Pump", "Hot Water", "Lighting", "Oven", "Others"))

analysisShortDT <- analysisShortDT[, c( "Total") := NULL]#We do not need this column

#Summarising the variation in demand
#We need to dcast it back into the original form
dt <- data.table::as.data.table(gather(analysisShortDT, eecaCircuit, totalEnergyWh, `Heat Pump`:`Others`, factor_key = TRUE))

dt <- dt[, totalEnergykWh := totalEnergyWh/1000]

myPlot <- ggplot2::ggplot(dt[year == 2015], aes((y=totalEnergykWh), x=eecaCircuit)) +
  geom_boxplot(aes(y=totalEnergykWh)) +
  labs(x = "Circuit", y= "Annual Energy consumption (kWh)") +
  #scale_x_continuous( limits = c(0,2000))+
  #facet_grid( season ~ year)+
  theme(text = element_text(family = "Cambria"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"))+
  theme(text = element_text(family = "Cambria")) 
myPlot

# kWh

# new DT ( do we need to?)
energyDT <- copy(powerDT)

#energyDT <- energyDT[, meanEnergyWh := meanPowerW/2]#Converting W into Wh energy 



#Conversion to energy Wh
energyDT <- energyDT[, meanEnergyWh := meanPowerW/2]
analysisDT <- copy(energyDT) # another new DT. Do we need to?
analysisDT <- analysisDT[, totalEnergyWh := sum(meanEnergyWh, na.rm = TRUE), 
                         keyby = .(linkID, season, eecaCircuit, year)]#Ignores NAs
analysisShortDT <- analysisDT[, .(totalEnergyWh = totalEnergyWh), 
                              keyby = .(linkID, season, eecaCircuit, year)] #For some reason by does not remove duplicates without crashing....

analysisShortDT <- unique(analysisShortDT)# Removing duplicates

# Calculating the contribution to annual energy consumption and turning off e-notation

options(scipen = 999)

analysisShortDT <- dcast(analysisShortDT, 
           linkID + year + season ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "totalEnergyWh") # take 1 household and turn data round

analysisShortDT <- analysisShortDT[, contribution_HW := `Hot water`/Total]
analysisShortDT <- analysisShortDT[, contribution_HP := `Heat Pump or Heating`/Total]
analysisShortDT <- analysisShortDT[, contribution_LI := `Lighting`/Total]
analysisShortDT <- analysisShortDT[, contribution_OV := `Oven`/Total]
analysisShortDT <- analysisShortDT[, contribution_XX := 1-((`Hot water` + `Heat Pump or Heating` + `Lighting` + `Oven`)/Total)]
# There are a lot of NA/Nan in the data. XXSensible to turn NaN into NA to receive more results??XX

analysisShortDT <- setnames(analysisShortDT, old=c("Heat Pump or Heating","Hot water", "Lighting", "Oven", "XX_Other"),
                        new=c("Heat Pump", "Hot Water", "Lighting", "Oven", "Others"))


#Summarising the variation in demand
#We need to dcast it back into the original form
dt <- data.table::as.data.table(gather(analysisShortDT, eecaCircuit, totalEnergyWh, `Heat Pump`:`Others`, factor_key = TRUE))

dt$season <- factor(dt$season, levels = c("Spring","Summer",
                                                    "Autumn", "Winter"))#This puts the seasons in the right order in the plot

Figure \@ref(fig:annualConsplot1) shows the distribution in seasonal energy consumption across all houses for Hot Water in 2015. Please note that Figure \@ref(fig:annualConsplot1) shows energy in kWh and not power demand in Watts. Results show that most energy consumption of Hot Water is in winter, followed by spring and autumn.

myPlot <- ggplot2::ggplot(dt[eecaCircuit =="Hot Water"], aes(x=totalEnergyWh/1000, colour=season)) +
  geom_density() +
  labs(x = "Seasonal energy (kWh)") +
  #scale_x_continuous( limits = c(0,2000))+
     geom_vline(aes(xintercept=mean(totalEnergyWh/1000)),
            color="blue", linetype="dashed", size=0.5)+
  #facet_grid( season ~ year)+

  theme(text = element_text(family = "Cambria"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"))
myPlot

In contrast, Figure \@ref(fig:annualConsplot2) presents the seasonal energy consumption for Heat Pumps and considers each dwelling and season. As expected, most energy consumption is in winter. In summer, total energy consumption is mostly 0 (appliance is turned off). Only a few observatons registered energy consumption in summer (cooling?).

myPlot <- ggplot2::ggplot(dt[eecaCircuit =="Heat Pump"], aes(x=totalEnergyWh/1000, colour=season)) +
  geom_density() +
  labs(x = "Seasonal energy (kWh)") +
  scale_x_continuous( limits = c(0,500))+
     geom_vline(aes(xintercept=mean(totalEnergyWh/1000)),
            color="blue", linetype="dashed", size=0.5)+
  #facet_grid( season ~ year)+

  theme(text = element_text(family = "Cambria"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"))
myPlot

Energy consumption of Lighting is particularly apparent in spring, autumn, and winter as presneted in Figure \@ref(fig:annualConsplot3). Less energy consumption is identified in summer.

myPlot <- ggplot2::ggplot(dt[eecaCircuit =="Lighting"], aes(x=totalEnergyWh/1000, colour=season)) +
  geom_density() +
  labs(x = "Seaonal energy (kWh)") +
  scale_x_continuous( limits = c(0,500))+
     geom_vline(aes(xintercept=mean(totalEnergyWh/1000)),
            color="blue", linetype="dashed", size=0.5)+
  #facet_grid( season ~ year)+

  theme(text = element_text(family = "Cambria"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"))
myPlot

For ovens, there is less seasonal variation in seasonal energy consumption as shown in Figure \@ref(fig:annualConsplot3). Furthermore, energy consumption in summer is less than in other seasons.

myPlot <- ggplot2::ggplot(dt[eecaCircuit =="Oven"], aes(x=totalEnergyWh/1000, colour=season)) +
  geom_density() +
  labs(x = "Seasonal energy (kWh)") +
  scale_x_continuous( limits = c(0,500))+
     geom_vline(aes(xintercept=mean(totalEnergyWh/1000)),
            color="blue", linetype="dashed", size=0.5)+
  #facet_grid( season ~ year)+

  theme(text = element_text(family = "Cambria"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"))
myPlot

Conservation Load Factor

#Part A

#The following calculates average power at 'peak' for each circuit as an average over all households  _XXKeep in mind that 11m data points are reduced to 6XX_
consloadDT <- copy(powerDT)#Utilise adjusted raw data

consloadDT <- consloadDT[season  == "Winter" &
             year == "2015" &
             obsHalfHour >= hms::as.hms("17:00:00") &
             obsHalfHour <= hms::as.hms("21:00:00")] #This is our definition of peak demand for this analysis regarding the conservation load factor

consloadDT <- consloadDT[, averageSeasonW := mean(meanPowerW, na.rm = TRUE), keyby =.(season, eecaCircuit, year, linkID)]#This gives us the average power at peak for each household and season and circuit

consloadDT <- consloadDT[, meanAverageSeasonW := mean(averageSeasonW, na.rm = TRUE), keyby =.(season, eecaCircuit, year)]#Average W over the households
consloadDT <- consloadDT[, .(meanAverageSeasonW =meanAverageSeasonW), keyby =.(season, eecaCircuit, year)]#Long to short dt
consloadDT<-  unique(consloadDT)# Removing duplicates XX_BE VERY CAREFUL WITH THIS_XX


#----------

#Part B


#The following calculates average energy consumption per year for each circuit as an average over all households _XXKeep in mind that 11m data points are reduced to 6XX_
EnergyConsDT <- copy(powerDT)#Utilise adjusted raw data to calculate average total energy consumption
EnergyConsDT <- EnergyConsDT[year == "2015"]#Only need to do this for 2015
EnergyConsDT <- EnergyConsDT[, EnergyWh := meanPowerW/2]#Creating energy column
EnergyConsDT <- EnergyConsDT[, totalEnergyWh := sum(EnergyWh, na.rm = TRUE), keyby =.(eecaCircuit, year, linkID)]#This gives us the total energy consumption for each circuit and household over all seasons
EnergyConsDT <- EnergyConsDT[, meantotalEnergyWh := mean(totalEnergyWh, na.rm = TRUE), keyby =.(eecaCircuit, year)]#Average Wh over the households
EnergyConsDT <- EnergyConsDT[, .(meantotalEnergyWh =meantotalEnergyWh), keyby =.(eecaCircuit, year)]#Long to short dt
EnergyConsDT<-  unique(EnergyConsDT)# Removing duplicates XX_BE VERY CAREFUL WITH THIS_XX

#----------

#Part C

#The following calculates the CLF based on the calculations above

consloadDT <- consloadDT[, meantotalEnergyWh := EnergyConsDT$meantotalEnergyWh]#Copy result in original conservation load table
consloadDT <- consloadDT[, consLoadFac := meanAverageSeasonW / (meantotalEnergyWh/1000)]#As defined by EECA







consloadCastDT <- dcast(consloadDT,
           season + year  ~ eecaCircuit, 
           fun = mean, # won't do anything as we should just have 1 value?
           value.var = "consLoadFac")
consloadCastDT<-consloadCastDT[, c( "Total") := NULL]#We do not need this column

consloadCastDT <- setnames(consloadCastDT, old=c("Heat Pump or Heating","Hot water", "Lighting", "Oven", "XX_Other"),
                        new=c("Heat Pump", "Hot Water", "Lighting", "Oven", "Others"))

consloadCastDT <- data.table::as.data.table(gather(consloadCastDT, eecaLabelContri, consLoadFac, `Heat Pump`:`Others`, factor_key = TRUE))#This brings our data into a nice form for the plot

Figure \@ref(fig:consLoadPlot) shows the Conservation Load Factor (CLF) for winter in 2015 as the devision of seasonal demand (W) and drawn energy (kWh) in this season as an average over all houses. The period of peak demand was set to 17:00-21:00 in winter. Based on this definition power was calculated as an average over all households. Results show that the CLF for Heat Pumps (0.56) is the highest followed by Lighting (0.45) and Oven (0.41). Please note that these numbers are averages over all households, dates, and times for 2015.

myPlot <- ggplot2::ggplot(consloadCastDT, aes(x= eecaLabelContri, fill = eecaLabelContri)) +
  geom_bar(aes(y = consLoadFac), stat = "identity")+
  theme(text = element_text(family = "Cambria"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(colour = "black"),
        legend.position = "none",
        axis.ticks.y = element_blank())+

    geom_text(aes(y = consLoadFac,
             label = sprintf("%0.2f",consLoadFac)),
            position = position_stack(vjust = .5))+

  #facet_grid(season ~ .) +
  #scale_y_continuous(breaks = c()) +
  #ggtitle("Contribution to peak demand (17:00-21:00) in 2015") +
  guides(fill=guide_legend(title="Circuit"))+
  labs(x='Circuit', y='Conservation Load Factor') 
myPlot

Summary

This report uses various definitions of peak demand to analyse the XXX

Data Annex

Figure \@ref(fig:winterPeakplot1) shows... for winter 2015 for the selected maximum load half-hours.

myPlot <- ggplot2::ggplot(resAllSeasonsDT[season == "Winter" & year == 2015], aes(x=Total)) +
  geom_histogram() +
  labs(x = "Mean total load per half hour") +
  scale_x_continuous( limits = c(0,13000))+
     geom_vline(aes(xintercept=mean(Total)),
            color="blue", linetype="dashed", size=0.5)+
  theme(text = element_text(family = "Cambria")) +
  ggtitle("Test") 
 # facet_grid(DOW ~ .) +
myPlot

Figure \@ref(fig:winterPeakplot2) shows...

myPlot <- ggplot2::ggplot(dataBucketDT[eecaCircuit  == "Heat Pump or Heating" & season == "Winter" & year == 2015],
                          aes(x=lubridate::date(r_dateTime_nz), group = lubridate::date(r_dateTime_nz))) +
  geom_boxplot(aes(y=maxPowerW, colour = maxPowerW), size=0.5, outlier.alpha = 0.1) +
  #scale_x_continuous( limits = c(0,13000))+
  # facet_grid(DOW ~ .) +
  theme(text = element_text(family = "Cambria")) 

myPlot

Figure \@ref(fig:winterPeakplot3) shows...

myPlot <- ggplot2::ggplot(resAllSeasonsDT[year==2015], aes(x=`Heat Pump or Heating`, colour = season)) +
  geom_density() +
  #scale_x_continuous( limits = c(0,2000000))+
     #geom_vline(aes(xintercept=mean(totalEnergyWh)),
           # color="blue", linetype="dashed", size=0.5)+
  theme(text = element_text(family = "Cambria")) +
  ggtitle("Test") 
myPlot

Figure \@ref(fig:winterPeakplot4) shows...

XX what is the gather for?

# put things where they are used

# dplyr so need to re-create as d.t
resDT <- data.table::as.data.table(gather(resAllSeasonsDT, eecaCircuit, maxPowerW, `Heat Pump or Heating`:`XX_Other`, factor_key = TRUE))

resDT <- resDT[, HP := NA]
resDT <- resDT[, `HP` := ifelse(`eecaCircuit`== "Heat Pump or Heating", maxPowerW, NA)]

myPlot <- ggplot2::ggplot(subset(resDT, year %in% c("2015") & season %in% c("Winter")),  fill = linkID) +
  geom_bar(aes(y= eecaCircuit, x = HP), stat = "identity", position = "dodge") +
  #geom_line(aes(y= V1)) +
  theme(text = element_text(family = "Cambria"))+
  #facet_grid(season ~ .) +
  ggtitle("XX") 
  #labs(x='Period', y='Cost in million NZD') 
myPlot

Categorised circuit labels

Table \@ref(tab:categoriseCircuitsTable) shows the categorised circuit labels used in the analysis.

kableExtra::kable(categoriseCircuitsTable, caption = "Categorised circuits - XX check carefully! Comment out when happy XX")  %>% # test result
  kable_styling()

Original data description

Descriptive statistics for aggregate half hourly power data for all households and all circuits:

skimr::skim(origPowerDT)

The following tables show descriptive statistics for the meanPowerW values for each circuit by household.

ids <- unique(origPowerDT$linkID)

for(hh in ids){
  # prints a lot of tables
    t <- data.table::cube(origPowerDT[linkID == hh], j = c(list(meanPowerW = mean(meanPowerW),
                                                              minPowerW = min(meanPowerW),
                                                              maxPowerW = max(meanPowerW),
                                                              nObs = .N)), by = c("circuit"))
    # NA in circuit column of ersults table = All
    t[, circuit := ifelse(is.na(circuit), "All", circuit)]
    print(kableExtra::kable(t, caption = paste0(hh, ": Mean of half-hourly mean power (W) by circuit type")))
}

Runtime

t <- proc.time() - startTime
elapsed <- t[[3]]

Analysis completed in r round(elapsed,2) seconds ( r round(elapsed/60,2) minutes) using knitr in RStudio with r R.version.string running on r R.version$platform.

R environment

R packages used

base R [@baseR]
bookdown [@bookdown]
data.table [@data.table]
ggplot2 [@ggplot2]
kableExtra [@kableExtra]
knitr [@knitr]
lubridate [@lubridate]
rmarkdown [@rmarkdown]

Session info

sessionInfo()

References

CfSOtago/GREENGridEECA documentation built on May 1, 2022, 8:08 p.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

CfSOtago/GREENGridEECA
Supports Analysis of NZ GREEN Grid data for EECA

In CfSOtago/GREENGridEECA: Supports Analysis of NZ GREEN Grid data for EECA

About

Citation

Report circulation:

License

History

Support

Introduction

Data

Peak Contribution

Winter (June - August) evenings 17:00 – 21:00

Region/network coincident peak

Hawke's Bay

Taranaki

Sample coincident peak

Annual Contribution

Conservation Load Factor

Summary

Data Annex

Categorised circuit labels

Original data description

Runtime

R environment

R packages used

Session info

References

R Package Documentation

Browse R Packages

We want your feedback!

CfSOtago/GREENGridEECA Supports Analysis of NZ GREEN Grid data for EECA

In CfSOtago/GREENGridEECA: Supports Analysis of NZ GREEN Grid data for EECA

About

Citation

Report circulation:

License

History

Support

Introduction

Data

Peak Contribution

Winter (June - August) evenings 17:00 – 21:00

Region/network coincident peak

Hawke's Bay

Taranaki

Sample coincident peak

Annual Contribution

Conservation Load Factor

Summary

Data Annex

Categorised circuit labels

Original data description

Runtime

R environment

R packages used

Session info

References

R Package Documentation

Browse R Packages

We want your feedback!

CfSOtago/GREENGridEECA
Supports Analysis of NZ GREEN Grid data for EECA