diff --git a/energy_demand/config_data/population_from_andrew.py b/energy_demand/config_data/population_from_andrew.py new file mode 100644 index 00000000..56704953 --- /dev/null +++ b/energy_demand/config_data/population_from_andrew.py @@ -0,0 +1,161 @@ +"""Download population projections from https://github.com/nismod/population/blob/master/README.md + + +Info +----- + +https://github.com/virgesmith/UKCensusAPI +https://www.nomisweb.co.uk/myaccount/webservice.asp +https://github.com/nismod/population +https://github.com/virgesmith/UKCensusAPI + +Steps +------ +1. optain nomis key +2. in cmd line: set NOMIS_API_KEY=XXX +3. ./setup.py install +4. run python script from command line + +https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationprojections/bulletins/nationalpopulationprojections/2015-10-29 + +Potential variants +------------------- +hhh: High population, +hpp: High fertility, +lll: Low population, +lpp: Low fertility, +php: High life expectancy, +pjp: Moderately high life expectancy, +pkp: Moderately low life expectancy, +plp: Low life expectancy, +pph: High migration, +ppl: Low migration, +ppp: Principal, +ppq: 0% future EU migration (non-ONS), +ppr: 50% future EU migration (non-ONS), +pps: 150% future EU migration (non-ONS), +ppz: Zero net migration + +Select for paper I +------------------- +Prinicpal projection ppp +Low mnigration ppl +High migration pph +""" +import pandas as pd +import population.nppdata as NPPData +import population.snppdata as SNPPData +import population.utils as utils + +reg_pop = True +extrapolate_pop = False +extrapolate_specific_scenario = True + +if extrapolate_specific_scenario: + ''' + Get extrapolated data for full time range for different ONS scenarios from 2015 - 2050 + # https://github.com/nismod/population/blob/master/doc/example_variant_ex.py + ''' + + # initialise the population modules + npp = NPPData.NPPData() + snpp = SNPPData.SNPPData() + + # 50 years, roughly half is extrapolated + years = range(2015, 2051) + + # start with an empty data frame + result_ppp = pd.DataFrame() + result_ppl = pd.DataFrame() + result_pph = pd.DataFrame() + + # loop over all the UK LAD (or LAD-equivalents) + for lad in snpp.data.GEOGRAPHY_CODE.unique(): + print("Collection population for region {}".format(lad)) + + # calculate the the variants + region_ppp = snpp.create_variant("ppp", npp, lad, years) + region_ppl = snpp.create_variant("ppl", npp, lad, years) + region_pph = snpp.create_variant("pph", npp, lad, years) + + # aggregate the calculated variants by age and gender + region_ppp = utils.aggregate(region_ppp, ["GENDER", "C_AGE"]) + region_ppl = utils.aggregate(region_ppl, ["GENDER", "C_AGE"]) + region_pph = utils.aggregate(region_pph, ["GENDER", "C_AGE"]) + + result_ppp.append(region_ppp, ignore_index=True) + result_ppl.append(region_ppl, ignore_index=True) + result_pph.append(region_pph, ignore_index=True) + + break + print("tt") + + # write out results + result_ppp.to_csv("C:/Users/cenv0553/mistral_population/__RESULTS/prinicpal_extrap_2050.csv", index=False) + result_ppl.to_csv("C:/Users/cenv0553/mistral_population/__RESULTS/lowmigration_extrap_2050.csv", index=False) + result_pph.to_csv("C:/Users/cenv0553/mistral_population/__RESULTS/highmigration_extrap_2050.csv", index=False) + print("FINISHED") +''' +if extrapolate_pop: + + # =================== + # Extrapolating SNPP for 2040 - 2050 + #https://github.com/nismod/population + # =================== + + # initialise the population modules + npp = NPPData.NPPData() + snpp = SNPPData.SNPPData() + + # get the first year where extrapolation is necessary + ex_start = snpp.max_year() + 1 + + # we extrapolate to 2050 + ex_end = 2050 + + # start with an empty data frame + result = pd.DataFrame() + + # loop over all the UK LAD (or LAD-equivalents) + for lad in snpp.data.GEOGRAPHY_CODE.unique(): + + # extrapolate + lad_ex = snpp.extrapolagg("GEOGRAPHY_CODE", npp, lad, range(ex_start, ex_end + 1)) + + # append to data + result = result.append(lad_ex, ignore_index=True) + + # write out results + result.to_csv("C:/Users/cenv0553/mistral_population/population/__RESULTS/snpp_extrap_2050.csv") + +if reg_pop: + + # initialise the population modules + npp = NPPData.NPPData() + snpp = SNPPData.SNPPData() + + # get the first year where extrapolation is necessary + ex_start = snpp.max_year() + 1 + + # we extrapolate to 2050 + ex_end = 2050 + + # start with an empty data frame + result = pd.DataFrame() + + # loop over all the UK LAD (or LAD-equivalents) + for lad in snpp.data.GEOGRAPHY_CODE.unique(): + + # get the total projected population for newcastle up to the SNPP horizon (2039) + lad_non_ex = snpp.aggregate("GEOGRAPHY_CODE", lad, range(2015, ex_start)) + + # extrapolate for another 25 years + lad_ex = snpp.extrapolagg("GEOGRAPHY_CODE", npp, lad, range(ex_end)) + + # append to data + result = result.append(lad_non_ex, ignore_index=True) + result = result.append(lad_ex, ignore_index=True) + + # write out results + result.to_csv("C:/Users/cenv0553/mistral_population/population/__RESULTS/snpp_extrap_2015_2050_ppp.csv") +''' \ No newline at end of file diff --git a/energy_demand/plotting/plotting_multiple_scenarios.py b/energy_demand/plotting/plotting_multiple_scenarios.py index fa87767b..311295df 100644 --- a/energy_demand/plotting/plotting_multiple_scenarios.py +++ b/energy_demand/plotting/plotting_multiple_scenarios.py @@ -637,7 +637,6 @@ def plot_tot_fueltype_y_over_time( # Set figure size fig = plt.figure(figsize=plotting_program.cm2inch(9, 8)) - ax = fig.add_subplot(1, 1, 1) y_scenario = {} @@ -845,9 +844,8 @@ def plot_tot_y_over_time( if plotshow: plt.show() - plt.close() - else: - plt.close() + + plt.close() def plot_radar_plots_average_peak_day( scenario_data, @@ -864,6 +862,8 @@ def plot_radar_plots_average_peak_day( name_spider_plot = os.path.join( fig_name, "spider_scenarios_{}.pdf".format(fueltype_to_model)) + fueltype_int = fueltypes[fueltype_to_model] + # ---------------- # Create base year peak load profile # Aggregate load profiles of all regions @@ -877,11 +877,12 @@ def plot_radar_plots_average_peak_day( print("-------Scenario: {} {}".format(scenario, fueltype_to_model)) base_yr = 2015 + # ------------------------ # Future year load profile + # ------------------------ all_regs_fueltypes_yh_by = np.sum(scenario_data[scenario]['results_every_year'][base_yr], axis=1) all_regs_fueltypes_yh_cy = np.sum(scenario_data[scenario]['results_every_year'][year_to_plot], axis=1) - fueltype_int = fueltypes[fueltype_to_model] # --------------------------- # Calculate load factors @@ -897,6 +898,30 @@ def plot_radar_plots_average_peak_day( scen_load_factor_fueltype_y_cy = load_factors.calc_lf_y(all_regs_fueltypes_yh_cy) load_factor_fueltype_y_cy.append(round(scen_load_factor_fueltype_y_cy[fueltype_int], fueltype_int)) + # ------------------------ + # Calculate share or space and water heating of total + # electrictiy demand in 2050 + # ------------------------ + enduses_to_agg = ['rs_space_heating', 'rs_water_heating', 'ss_space_heating', 'ss_water_heating', 'is_space_heating'] + + aggregated_enduse_fueltype = np.zeros((365, 24)) + aggregated_enduse_fueltype_by = np.zeros((365, 24)) + total_demand = np.zeros((365, 24)) + for enduse in scenario_data[scenario]['results_enduse_every_year'][2050].keys(): + if enduse in enduses_to_agg: + aggregated_enduse_fueltype += scenario_data[scenario]['results_enduse_every_year'][2050][enduse][fueltype_int] + aggregated_enduse_fueltype_by += scenario_data[scenario]['results_enduse_every_year'][2015][enduse][fueltype_int] + total_demand += scenario_data[scenario]['results_enduse_every_year'][2050][enduse][fueltype_int] + else: + total_demand += scenario_data[scenario]['results_enduse_every_year'][2050][enduse][fueltype_int] + + # Total demand + selected_enduses_total_peak_day = round(np.max(aggregated_enduse_fueltype[peak_day_nr_cy]), 1) + all_enduses_total_peak_day = round(np.max(total_demand[peak_day_nr_cy]), 1) + selected_enduses_peak_by = round(np.max(aggregated_enduse_fueltype_by), 1) + print("p_all_enduses_total_peak_day " + str(all_enduses_total_peak_day)) + print("p_selected_enduses_total_peak_day " + str(selected_enduses_total_peak_day)) + # ---------- # Calculate change in peak # ---------- @@ -905,11 +930,14 @@ def plot_radar_plots_average_peak_day( diff_max_h = round(((100 / by_max_h) * cy_max_h) - 100, 1) - label_max_h = " by: {} cy: {} d: {}, scen: {}".format( + label_max_h = " by: {} cy: {} d: {}, scen: {} ed_peakh_heating: {} ed_tot_peah_ {} by_heating: {}".format( round(by_max_h, 2), round(cy_max_h, 1), round(diff_max_h, 1), - scenario,) + scenario, + selected_enduses_total_peak_day, + all_enduses_total_peak_day, + selected_enduses_peak_by) list_diff_max_h.append(label_max_h) diff --git a/energy_demand/plotting/plotting_results.py b/energy_demand/plotting/plotting_results.py index b9c863c9..93c39906 100644 --- a/energy_demand/plotting/plotting_results.py +++ b/energy_demand/plotting/plotting_results.py @@ -2219,34 +2219,26 @@ def plot_radar_plot_multiple_lines( color_line, alpha=0.05) - font_additional_info = plotting_styles.font_info(size=5) + # ------------ + # Write outs + # ------------ + font_additional_info = plotting_styles.font_info(size=4) for cnt, entry in enumerate(list_diff_max_h): plt.text( - 0.25, + 0.15, 0 + cnt/50, entry, fontdict=font_additional_info, transform=plt.gcf().transFigure) - # ------------ - # Title - # ------------ - font_additional_info = plotting_styles.font_info(size=4) - #plt.title( - # title_info, - # loc='left', - # fontdict=font_additional_info) - # ------------ # Legend # ------------ plt.legend( ncol=1, - bbox_to_anchor=(0.5, -0.1), - prop={ - 'family': 'arial', - 'size': 4}, + bbox_to_anchor=(0.2, -0.1), + prop={'size': 4}, frameon=False) plt.savefig(fig_name) @@ -2267,8 +2259,8 @@ def plt_one_fueltype_multiple_regions_peak_h( Arguments --------- """ - # Set figure size - fig = plt.figure(figsize=plotting_program.cm2inch(14, 8)) + fig = plt.figure( + figsize=plotting_program.cm2inch(14, 8)) ax = fig.add_subplot(1, 1, 1) diff --git a/energy_demand/plotting/result_mapping.py b/energy_demand/plotting/result_mapping.py index 76267b25..f010dda4 100644 --- a/energy_demand/plotting/result_mapping.py +++ b/energy_demand/plotting/result_mapping.py @@ -15,92 +15,6 @@ from energy_demand.basic import basic_functions from energy_demand.technologies import tech_related -# Old reasonable function to generate symetric values -'''def get_reasonable_bin_values( - data_to_plot, - increments=10 - ): - """Get reasonable bin values - """ - def round_down(num, divisor): - return num - (num%divisor) - - max_val = max(data_to_plot) - min_val = min(data_to_plot) - - # Positive values - if min_val >= 0: - - # Calculate number of classes - classes = max_val / increments - - # Round down and add one class for larger values - nr_classes = round_down(classes, increments) + 1 - - if nr_classes > 9: - raise Exception("Nr of classes is too big") - - # Classes - min_class = round_down(min_val, increments) # Minimum class - max_class = round_down(max_val, increments) + increments # Maximum class - - # Bin with classes - bins = range(min_class, max_class, increments) - else: - #lager negative values - - # must be of uneven length containing zero - largest_min = abs(min_val) - largest_max = abs(max_val) - - nr_min_class = round_down(abs(min_val), increments) / increments - - if max_val < 0: - nr_pos_classes = nr_min_class - else: - nr_pos_classes = round_down(max_val, increments) / increments - - if nr_min_class > nr_pos_classes: - # Negative has more classes - nr_of_classes_symetric = nr_min_class - else: - nr_of_classes_symetric = nr_pos_classes - - # Number of classes - symetric_value = nr_of_classes_symetric * increments - min_class_value = int(symetric_value * -1) - max_class_value = int(symetric_value) - - # Negative classes - if min_class_value / increments == 0: - neg_classes = [increments * -1] - elif min_class_value / increments == 1: - neg_classes = [increments * -1] - else: - neg_classes = list(range(min_class_value, 0, increments)) - - if max_class_value / increments == 0: - pos_classes = [increments] - elif max_class_value / increments == 1: - pos_classes = [increments] #only one class - else: - pos_classes = list(range(increments, max_class_value + increments, increments)) - - bins = neg_classes + [0] + pos_classes - - logging.info("______________________________________") - logging.info("largest_min " + str(largest_min)) - logging.info("largest_max " + str(largest_max)) - logging.info(min_val) - logging.info(max_val) - logging.info("symetric_value " + str(symetric_value)) - logging.info("min_class_value " + str(min_class_value)) - logging.info("max_class_value " + str(max_class_value)) - logging.info("bins " + str(bins)) - - return bins -''' - def get_reasonable_bin_values( data_to_plot, increments=10 @@ -133,8 +47,8 @@ def round_down(num, divisor): min_class = increments # Bin with classes - logging.info("vv {} {}".format(max_val, min_val)) - logging.info("pos {} {} {}".format(min_class, max_class, increments)) + #logging.info("vv {} {}".format(max_val, min_val)) + #logging.info("pos {} {} {}".format(min_class, max_class, increments)) bins = list(range(int(min_class), int(max_class), int(increments))) else: logging.info("Neg") @@ -145,18 +59,12 @@ def round_down(num, divisor): largest_max = abs(max_val) nr_min_class = round_down(abs(min_val), increments) / increments - + if max_val < 0: nr_pos_classes = 0 else: nr_pos_classes = round_down(max_val, increments) / increments + 1 - '''if nr_min_class > nr_pos_classes: - # Negative has more classes - nr_of_classes_symetric = nr_min_class - else: - nr_of_classes_symetric = nr_pos_classes''' - # Number of classes #symetric_value = nr_of_classes_symetric * increments min_class_value = int(nr_min_class * -1) * increments @@ -170,7 +78,7 @@ def round_down(num, divisor): neg_classes = [increments * -1] else: neg_classes = list(range(min_class_value, 0, increments)) - + if max_class_value == 0: pos_classes = [] elif max_class_value / increments == 0: @@ -179,9 +87,15 @@ def round_down(num, divisor): pos_classes = [increments] #only one class else: pos_classes = list(range(increments, max_class_value + increments, increments)) - + bins = neg_classes + [0] + pos_classes + # --- + # Test that maximum 9 classes + # --- + if len(bins) > 8: + raise Exception("Too many bin classes defined") + return bins def user_defined_classification( @@ -400,7 +314,6 @@ def re_classification( """ # Create the list of bin labels and the # list of colors corresponding to each bin - logging.info("EE " + str(bins)) if max(bins) <=0: nr_of_classes = int(len(bins)) + 1 #only negative classes elif min(bins) >= 0: @@ -409,23 +322,20 @@ def re_classification( nr_of_classes = int(len(bins) -1) + 2 #class on top and bottom and get rid of 0 bin_labels = [] - - logging.info("nr_of_classes" + str(nr_of_classes)) - #for idx, _ in enumerate(bins): - # val_bin = idx / (len(bins) - 1.0) - # bin_labels.append(val_bin) + for idx in range(nr_of_classes): val_bin = idx / (nr_of_classes - 1.0) bin_labels.append(val_bin) + # ---------------------- # Add white bin and color # ---------------------- color_list_copy = copy.copy(color_list) bin_labels_copy = copy.copy(bin_labels) - logging.info("KUH " + str(color_zero)) - logging.info(color_list_copy) - logging.info(bin_labels_copy) + #logging.info("KUH " + str(color_zero)) + #logging.info(color_list_copy) + #logging.info(bin_labels_copy) # Add zero color in color list if a min_plus map if color_zero != False: insert_pos = 1 @@ -435,17 +345,17 @@ def re_classification( bin_labels_copy.insert(insert_pos, float(placeholder_zero_color)) else: pass - logging.info("KUH 2") - logging.info(color_list_copy) - logging.info(bin_labels_copy) + #logging.info("KUH 2") + #logging.info(color_list_copy) + #logging.info(bin_labels_copy) # Create the custom color map color_bin_match_list = [] for lbl, color in zip(bin_labels_copy, color_list_copy): color_bin_match_list.append((lbl, color)) - logging.info("color_bin_match_list: " + str(color_bin_match_list)) - logging.info(bins) - logging.info(bin_labels_copy) + ##logging.info("color_bin_match_list: " + str(color_bin_match_list)) + ##logging.info(bins) + ##logging.info(bin_labels_copy) if 0 in bins: pass @@ -453,17 +363,16 @@ def re_classification( logging.info("TT") #bins.insert(0, 0) #TODO else: - logging.info("aa") bins.insert(len(bins), 0) # Add zero at the end - logging.info("Newbin") - logging.info(bins) + #logging.info("Newbin") + #logging.info(bins) cmap = LinearSegmentedColormap.from_list( 'mycmap', color_bin_match_list) - logging.info("cmap " + str(cmap)) + #logging.info("cmap " + str(cmap)) # Reclassify lad_geopanda_shp['reclassified'] = lad_geopanda_shp[field_to_plot].apply( @@ -1082,11 +991,10 @@ def create_geopanda_files( # =============================================== # Differences in percent per enduse and year (y) # =============================================== - if plot_crit_dict['plot_differences_p'] and year == 2050: #base_yr: #TODO TODO + if plot_crit_dict['plot_differences_p'] and year > 2015: field_name = 'y_diff_p_{}-{}_{}'.format( base_yr, year, fueltype_str) - logging.info("===========OKO " + str(field_name)) - print("OKO " + str(field_name)) + #logging.info("===========field_name " + str(field_name)) # Calculate yearly sums yearly_sum_gwh_by = np.sum( @@ -1122,29 +1030,27 @@ def create_geopanda_files( if nan_entry: continue - bins_increments = 15 + # ---- + # CAlculate classes for manual classification + # ---- + bins_increments = 10 + bins = get_reasonable_bin_values( data_to_plot=list(data_to_plot.values()), increments=bins_increments) - - logging.info("Min {} Max {}".format( - min(list(data_to_plot.values())), - max(list(data_to_plot.values())))) - logging.info("FIRST BINS: " + str(bins)) - #bins = [-4, -2, 0, 2, 4] +# must be of uneven length containing zero + + #logging.info("Min {} Max {}".format( + # min(list(data_to_plot.values())), + # max(list(data_to_plot.values())))) + #logging.info("FIRST BINS: " + str(bins)) + #bins = [-4, -2, 0, 2, 4] #bins = [-40, -30, -20, -10, 0, 10, 20, 30, 40] color_list, color_prop, user_classification, color_zero = colors_plus_minus_map( bins=bins, color_prop='qualitative', color_order=True) - - logging.info("KURT") - logging.info(color_list) - logging.info(color_prop) - logging.info(user_classification) - logging.info("_----------------") - user_classification = True #NEW TODO + #user_classification = True #NEW TODO # Plot difference in % per fueltype of total fuel (y) plot_lad_national( diff --git a/energy_demand/processing/multiple_scenarios.py b/energy_demand/processing/multiple_scenarios.py index 3cd9a46c..ab10cb0e 100644 --- a/energy_demand/processing/multiple_scenarios.py +++ b/energy_demand/processing/multiple_scenarios.py @@ -31,8 +31,8 @@ def process_result_multi_scen(path_to_scenarios, path_shapefile_input): # Execute rusult processing for every scenario #process_result_multi_scen("C:/Users/cenv0553/ed/results/scen_A", os.path.abspath('C:/Users/cenv0553/ED/data/_raw_data/C_LAD_geography/same_as_pop_scenario/lad_2016_uk_simplified.shp')) #process_result_multi_scen("C:/Users/cenv0553/ed/results/scen_B", os.path.abspath('C:/Users/cenv0553/ED/data/_raw_data/C_LAD_geography/same_as_pop_scenario/lad_2016_uk_simplified.shp')) -#process_result_multi_scen("C:/Users/cenv0553/ed/results/scen_C", os.path.abspath('C:/Users/cenv0553/ED/data/_raw_data/C_LAD_geography/same_as_pop_scenario/lad_2016_uk_simplified.shp')) -process_result_multi_scen("C:/Users/cenv0553/ed/results/scen_D", os.path.abspath('C:/Users/cenv0553/ED/data/_raw_data/C_LAD_geography/same_as_pop_scenario/lad_2016_uk_simplified.shp')) +process_result_multi_scen("C:/Users/cenv0553/ed/results/scen_C", os.path.abspath('C:/Users/cenv0553/ED/data/_raw_data/C_LAD_geography/same_as_pop_scenario/lad_2016_uk_simplified.shp')) +#process_result_multi_scen("C:/Users/cenv0553/ed/results/scen_D", os.path.abspath('C:/Users/cenv0553/ED/data/_raw_data/C_LAD_geography/same_as_pop_scenario/lad_2016_uk_simplified.shp')) #process_result_multi_scen("C:/Users/cenv0553/ed/results/Fig_11", os.path.abspath('C:/Users/cenv0553/ED/data/_raw_data/C_LAD_geography/same_as_pop_scenario/lad_2016_uk_simplified.shp')) diff --git a/energy_demand/processing/scenario_charts.py b/energy_demand/processing/scenario_charts.py index 9e6cb030..05bdf881 100644 --- a/energy_demand/processing/scenario_charts.py +++ b/energy_demand/processing/scenario_charts.py @@ -187,6 +187,6 @@ def process_scenarios(path_to_scenarios, year_to_model=2015): process_scenarios(os.path.abspath("C:/Users/cenv0553/ed/results/scen_C")) process_scenarios(os.path.abspath("C:/Users/cenv0553/ed/results/scen_D"))''' -#process_scenarios(os.path.abspath("C:/Users/cenv0553/ed/results/Fig_08_09")) +process_scenarios(os.path.abspath("C:/Users/cenv0553/ed/results/Fig_08_09")) #process_scenarios(os.path.abspath("C:/Users/cenv0553/ed/results/Fig_11")) -process_scenarios(os.path.abspath("C:/Users/cenv0553/ed/results/Fig_13")) # Color scheme Fig 13 +#process_scenarios(os.path.abspath("C:/Users/cenv0553/ed/results/Fig_13")) # Color scheme Fig 13