Skip to content

Commit

Permalink
Add segment HDH/CDH for grid region and Update to use time period pop…
Browse files Browse the repository at this point in the history 
…ulation
  • Loading branch information
@mengqi-z
mengqi-z committed Mar 7, 2024
1 parent 9cfc798 commit 7c2b2ac
Show file tree
Hide file tree
Showing 4 changed files with 240 additions and 36 deletions.
2 changes: 1 addition & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ vignettes/*.pdf

# R Environment Variables
.Renviron
output
output*
inst/doc
docs/
tests/testthat/output/
Expand Down
103 changes: 95 additions & 8 deletions R/diagnostics.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,27 @@ diagnostics <- function(hdcd_segment = tibble::tibble(),
'Nov_day', 'Nov_night', 'Dec_day', 'Dec_night',
'superpeak')

if(any(grepl('degree-hours', unique(hdcd_segment$unit)))){

# calculate the number of hours in each segment by state
segment_hours <- helios::segment_map_utc %>%
dplyr::group_by(subRegion, segment) %>%
dplyr::summarise(segment_hours = n()) %>%
dplyr::ungroup()

if(any(grepl('grid', unique(hdcd_segment$subRegion)))) {
# calculate the number of hours in each segment by grid region
segment_hours <- segment_hours %>%
dplyr::left_join(helios::mapping_states_gridregion, by = 'subRegion') %>%
dplyr::select(-subRegion) %>%
dplyr::rename(subRegion = grid_region) %>%
dplyr::distinct()
}


}


if(any(unique(hdcd_segment$segment) %in% segment_levels)) {

# Check if the outputs cover a full year
Expand Down Expand Up @@ -89,6 +110,7 @@ diagnostics <- function(hdcd_segment = tibble::tibble(),
# Individual Years
for(year_i in (hdcd_comb_diagnostics$year) %>% unique()) {

# plot heating and cooling degree-hours ------------------------------
data_plot <- hdcd_comb_diagnostics %>%
dplyr::filter(year == year_i) %>%
dplyr::mutate(segment = factor(segment, levels = segment_levels))
Expand Down Expand Up @@ -119,10 +141,46 @@ diagnostics <- function(hdcd_segment = tibble::tibble(),
height = 15)

print(paste0('Diagnostic figure saved as : ', filename_diagnostics_i))

# plot heating and cooling degree-hours **load** ---------------------
data_plot <- hdcd_comb_diagnostics %>%
dplyr::filter(year == year_i) %>%
dplyr::left_join(segment_hours, by = c('subRegion', 'segment')) %>%
dplyr::mutate(value = value / segment_hours,
segment = factor(segment, levels = segment_levels))

# plot
p <- ggplot2::ggplot(data = data_plot) +
ggplot2::aes(x = segment, y = value, group = heatcool) +
ggplot2::geom_line(ggplot2::aes(color = heatcool)) +
ggplot2::facet_wrap(subRegion ~ ., scales = 'free_y') +
ggplot2::ggtitle(paste0('HDCD Load at GCAM-USA Dispatch Segment in ', year_i)) +
ggplot2::ylab('Degree Load') +
ggplot2::scale_color_manual(values = c('heat' = '#1AB2FF',
'cool' = '#E61A33')) +
ggplot2::scale_x_discrete(drop = FALSE) +
ggplot2::theme_bw() +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90,
vjust = 0.5))

# create plot name
filename_diagnostics_i <- file.path(
folder_diagnostics,
helios::create_name(c('segment_load', year_i, name_append), 'png'))

# save plot
ggplot2::ggsave(p,
filename = filename_diagnostics_i,
width = 25,
height = 15)

print(paste0('Diagnostic figure saved as : ', filename_diagnostics_i))
}

# combined years with color gradients with free scale
if(length(unique(hdcd_comb_diagnostics$year)) > 1){

# Plot all years for HDCD degree-hours by segment --------------------
data_plot <- hdcd_comb_diagnostics %>%
dplyr::mutate(segment = factor(segment, levels = segment_levels))

Expand Down Expand Up @@ -165,21 +223,50 @@ diagnostics <- function(hdcd_segment = tibble::tibble(),

}

}

# combined years with color gradients with fixed scale
if(length(unique(hdcd_comb_diagnostics$year)) > 1){
# Plot all years for HDCD degree **load** by segment -----------------
data_plot <- hdcd_comb_diagnostics %>%
dplyr::mutate(segment = factor(segment, levels = segment_levels))
dplyr::left_join(segment_hours, by = c('subRegion', 'segment')) %>%
dplyr::mutate(value = value / segment_hours,
segment = factor(segment, levels = segment_levels))

n_color <- length(unique(hdcd_comb_diagnostics$year))
pal_hd <- colorRampPalette(RColorBrewer::brewer.pal(9, 'YlOrRd'))
pal_cd <- colorRampPalette(RColorBrewer::brewer.pal(9, 'YlGnBu'))
pal <- c(rev(pal_hd(n_color)), rev(pal_cd(n_color)))
for (scale_i in c('free_y', 'fixed')){

scale_name <- dplyr::case_when(scale_i == 'free_y' ~ 'freeScale',
scale_i == 'fixed' ~ 'fixedScale')

# plot
p <- ggplot2::ggplot(data = data_plot) +
ggplot2::geom_line(ggplot2::aes(x = segment, y = value,
group = interaction(year, heatcool),
color = interaction(year, heatcool))) +
ggplot2::facet_wrap(subRegion ~ ., scales = scale_i) +
ggplot2::ggtitle(paste0('HDCD Load at GCAM-USA Dispatch Segment for ', hdcd_comb_year_range)) +
ggplot2::ylab('Degree Load') +
ggplot2::scale_color_manual(values = pal,
guide = ggplot2::guide_legend(title = 'HDCD (All Years)')) +
ggplot2::scale_x_discrete(drop = FALSE) +
ggplot2::theme_bw() +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90,
vjust = 0.5))

# create plot name
filename_diagnostics_i <- file.path(
folder_diagnostics,
helios::create_name(c('segment_load_allYears', scale_name, name_append), 'png'))

# save plot
ggplot2::ggsave(p,
filename = filename_diagnostics_i,
width = 25,
height = 15)
print(paste0('Diagnostic figure saved as : ', filename_diagnostics_i))

}

}


} else {
message(paste0('Data is less than ', min_diagnostic_months, ' months. No diagnostic plots by dispatch segment are created.'))
} # end of if(length(unique(hdcd_segment$segment)) >= min_diagnostic_months * 2)
Expand Down
159 changes: 132 additions & 27 deletions R/hdcd.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ hdcd <- function(ncdf = NULL,

# for IM3 grid region analysis
if(im3_analysis){
hdcd_comb_gridregion_gcam <- tibble::tibble()
hdcd_comb_gridregion_monthly <- tibble::tibble()
hdcd_comb_gridregion_annual <- tibble::tibble()
}
Expand Down Expand Up @@ -180,11 +181,48 @@ hdcd <- function(ncdf = NULL,

population_j = population[[j]]

# read population based on the population data type
# output from wrf resolution: ['ID', 'subRegion', 'lat', 'lon', 'year', 'value' ]
# output from 1/8th degree pop data: ['ID', 'subRegion', 'lat', 'lon', 'year', 'value' ]
population_j_grid <- helios::read_population(file = population_j,
time_periods = time_periods)
if (im3_analysis) {

# use the population at each time period for the HDCD years within that period
# (1) year 2020 for period 2020
# (2) year 2021 - 2025 for period 2025, ..., year 2091-2095 for period 2095
# (3) year 2096 - 2099 for period 2100
# get the trend-representative year for each 5-year period
periods <- c(as.Date("2020-01-01"),
as.Date("2021-01-01"),
seq(as.Date("2026-01-01"), as.Date("2101-01-01"), by= "5 years"))

mapping_time_period <- data.frame(year = time_periods) %>%
dplyr::mutate(date = as.Date(paste0(year, '-01-01')),
period = cut(date, breaks = periods, right = FALSE, labels = periods[2:length(periods)]),
period = lubridate::year(as.Date(period) - 1)) %>%
dplyr::select(-date)

# read population based on the population data type
# output from wrf resolution: ['ID', 'subRegion', 'lat', 'lon', 'year', 'value' ]
# output from 1/8th degree pop data: ['ID', 'subRegion', 'lat', 'lon', 'year', 'value' ]
population_j_grid <- helios::read_population(file = population_j,
time_periods = unique(mapping_time_period$period))

# expand population data from 5-year period to annual based on the mapping_time_period
# the annual population for years within the period will be the same as the population for the period
# E.g., population from years 2021 - 2025 uses population in period 2025
population_j_grid <- population_j_grid %>%
tidyr::expand(tidyr::nesting(lat, lon),
year = time_periods) %>%
dplyr::left_join(mapping_time_period, by = 'year') %>%
dplyr::left_join(population_j_grid, by = c('lat', 'lon', 'period' = 'year')) %>%
dplyr::select(-period)

} else {

# read population based on the population data type
# output from wrf resolution: ['ID', 'subRegion', 'lat', 'lon', 'year', 'value' ]
# output from 1/8th degree pop data: ['ID', 'subRegion', 'lat', 'lon', 'year', 'value' ]
population_j_grid <- helios::read_population(file = population_j,
time_periods = time_periods)

}

population_j_grid <- helios::match_grids(from_df = population_j_grid,
to_df = ncdf_pivot,
Expand All @@ -204,8 +242,6 @@ hdcd <- function(ncdf = NULL,
spatial = spatial)




# Create population weighted raster if any population years in ncdf years
if (any(unique(population_j_grid$year) %in% years)) {

Expand Down Expand Up @@ -362,26 +398,33 @@ hdcd <- function(ncdf = NULL,
#......................
# Step 6a: Aggregate over Segments, monthly and annual
#......................
#
# temporal_subset <- data.frame(ncdf_times = ncdf_times[index_subset],
# x = paste0('X', index_subset)) %>%
# dplyr::mutate(datetime = as.POSIXct(ncdf_times, format = '%Y-%m-%d_%H:%M:%S', tz = 'UTC')) %>%
# dplyr::mutate(year = lubridate::year(datetime),
# month = lubridate::month(datetime),
# day = lubridate::day(datetime),
# hour = lubridate::hour(datetime),
# timezone = lubridate::tz(datetime)) %>%
# dplyr::mutate(month = dplyr::if_else(month < 10, paste0('0', month), paste0(month)),
# day = dplyr::if_else(day < 10, paste0('0', day), paste0(day)),
# hour = dplyr::if_else(hour < 10, paste0('0', hour), paste0(hour)) )
#
# # hdcd segments
# hdcd_region_segments <- hdcd_region %>%
# dplyr::left_join(temporal_subset, by = c('datetime', 'year')) %>%
# dplyr::left_join(helios::segment_map_utc,
# by = c('subRegion', 'month', 'day', 'hour')) %>%
# dplyr::group_by(region, subRegion, year, segment, HDCD) %>%
# dplyr::summarise(value = sum(value, na.rm = T))

temporal_subset_gridregion <- data.frame(ncdf_times = ncdf_times[index_subset],
x = paste0('X', index_subset)) %>%
dplyr::mutate(datetime = as.POSIXct(ncdf_times, format = '%Y-%m-%d_%H:%M:%S', tz = 'UTC')) %>%
dplyr::mutate(year = lubridate::year(datetime),
month = lubridate::month(datetime),
day = lubridate::day(datetime),
hour = lubridate::hour(datetime),
timezone = lubridate::tz(datetime)) %>%
dplyr::mutate(month = dplyr::if_else(month < 10, paste0('0', month), paste0(month)),
day = dplyr::if_else(day < 10, paste0('0', day), paste0(day)),
hour = dplyr::if_else(hour < 10, paste0('0', hour), paste0(hour)) )

# segment map for grid regions
segment_map_utc_gridregion <- helios::segment_map_utc %>%
dplyr::left_join(helios::mapping_states_gridregion, by = 'subRegion') %>%
dplyr::select(-subRegion) %>%
dplyr::distinct() %>%
dplyr::rename(subRegion = grid_region)

# hdcd segments **degree hours**
hdcd_gridregion_segments <- hdcd_gridregion %>%
dplyr::left_join(temporal_subset, by = c('datetime', 'year')) %>%
dplyr::left_join(segment_map_utc_gridregion,
by = c('subRegion', 'month', 'day', 'hour')) %>%
dplyr::group_by(region, subRegion, year, segment, HDCD) %>%
dplyr::summarise(value = sum(value, na.rm = T))

# hdcd monthly **degree hours**
hdcd_gridregion_monthly <- hdcd_gridregion %>%
Expand Down Expand Up @@ -427,6 +470,29 @@ hdcd <- function(ncdf = NULL,
!((value > 0) & grepl('heat', thermal.building.service.input, ignore.case = T))) %>%
dplyr::select(-HDCD)


# for IM3 analysis
if(im3_analysis){

# create building service structure for grid region
L2441.HDDCDD_Fixed_gcamusa_seg_gridregion <- helios::L2441.HDDCDD_Fixed_gcamusa_seg %>%
dplyr::left_join(helios::mapping_states_gridregion, by = 'subRegion') %>%
dplyr::select(-subRegion) %>%
dplyr::rename(subRegion = grid_region) %>%
dplyr::distinct()

# merge with segment table
hdcd_gridregion_bld <- hdcd_gridregion_segments %>%
dplyr::left_join(L2441.HDDCDD_Fixed_gcamusa_seg_gridregion,
by = c('subRegion', 'segment')) %>%
# Remove columns with -ve hdcd/cooling and +ve/heating
dplyr::filter(
!((value < 0) & grepl('cool', thermal.building.service.input, ignore.case = T)),
!((value > 0) & grepl('heat', thermal.building.service.input, ignore.case = T))) %>%
dplyr::select(-HDCD)

}

} else {
# only gcam_us49 can have dispatch segment because
# we don't have day, night segement for other regions
Expand Down Expand Up @@ -578,6 +644,24 @@ hdcd <- function(ncdf = NULL,

hdcd_name <- 'hdhcdh'

# for im3 analysis
if(im3_analysis){

hdcd_comb_gridregion_gcam <- hdcd_comb_gridregion_gcam %>%
dplyr::bind_rows(hdcd_gridregion_bld %>%
dplyr::mutate(unit = 'Fahrenheit degree-hours')) %>%
dplyr::ungroup() %>%
dplyr::group_by(region, subRegion, year, segment, gcam.consumer,
nodeInput, building.node.input,
thermal.building.service.input, unit) %>%
dplyr::summarize(value = sum(value, na.rm = T)) %>%
dplyr::ungroup() %>%
dplyr::filter(!is.na(subRegion),
!is.na(year),
!is.na(segment))

}

} else if ((model_timestep == 'daily' & dispatch_segment == FALSE) |
(model_timestep == 'hourly' & dispatch_segment == FALSE)){
hdcd_comb_gcam <- hdcd_comb_gcam %>%
Expand Down Expand Up @@ -614,6 +698,26 @@ hdcd <- function(ncdf = NULL,
}
}

# for IM3 analysis
if(im3_analysis){

if(i == length(ncdf)){

year_min_i <- min(hdcd_comb_gridregion_gcam$year, na.rm = T)
year_max_i <- max(hdcd_comb_gridregion_gcam$year, na.rm = T)

filename_i_gridregion <- file.path(
folder,
helios::create_name(c('hdhcdh', model, year_min_i, year_max_i, 'gcam', 'gridregion', name_append), 'csv'))

if(save){
data.table::fwrite(hdcd_comb_gridregion_gcam, file = filename_i_gridregion)
print(paste0('File saved as : ', filename_i_gridregion))
}
}

}

}

#......................
Expand Down Expand Up @@ -791,6 +895,7 @@ hdcd <- function(ncdf = NULL,
invisible(list(hdcd_comb_gcam = hdcd_comb_gcam,
hdcd_comb_monthly = hdcd_comb_monthly,
hdcd_comb_annual = hdcd_comb_annual,
hdcd_comb_gridregion_gcam = hdcd_comb_gridregion_gcam,
hdcd_comb_gridregion_monthly = hdcd_comb_gridregion_monthly,
hdcd_comb_gridregion_annual = hdcd_comb_gridregion_annual))

Expand Down
12 changes: 12 additions & 0 deletions inst/extras/dev_tests.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,3 +240,15 @@ helios::diagnostics(
folder = file.path(getwd(), 'output'),
name_append = 'combined_hdhcdh_monthly_gridregion'
)

# Test diagnostics with NERSC output on segment HDHCDH at US49
hdhcdh_us49 <- data.table::fread(
'C:/WorkSpace/IM3/helios/hddcdd/nersc/combined_outputs_hdcd_rcp45cooler_ssp3/combined_hdhcdh_2020_2099_gcam_us49.csv'
)

helios::diagnostics(
hdcd_segment = hdhcdh_us49,
min_diagnostic_months = 6,
folder = file.path(getwd(), 'output'),
name_append = 'combined_hdhcdh_gcam_us49'
)

0 comments on commit 7c2b2ac

Please sign in to comment.