From b3bb4c62d0457ef5621868f91a94db2d5a90a00a Mon Sep 17 00:00:00 2001 From: Sven Eggimann Date: Fri, 28 Sep 2018 11:53:01 +0100 Subject: [PATCH] fixed bug in line smoothing --- energy_demand/main.py | 2 +- energy_demand/plotting/_scrap.py | 38 +---- .../plotting/basic_plot_functions.py | 80 +++++++--- .../plotting/fig_load_profile_dh_multiple.py | 25 ++-- .../plotting/fig_total_demand_peak.py | 139 ++++++++++++++---- energy_demand/plotting/result_mapping.py | 1 + .../result_processing/scenario_charts.py | 2 - .../result_processing/single_scenario.py | 65 ++++---- 8 files changed, 222 insertions(+), 130 deletions(-) diff --git a/energy_demand/main.py b/energy_demand/main.py index 52bd9533..51de5016 100644 --- a/energy_demand/main.py +++ b/energy_demand/main.py @@ -10,7 +10,7 @@ # MAKE SIMLPLE TABLE FOR READING IN FUELS # correction factors on LAD level disaggregation? (load non-residential demand) # Improve plotting and processing (e.g. saisonal plots) - +# test if smoothing and plots # Weather station cleaning: Replace days with missing values Note diff --git a/energy_demand/plotting/_scrap.py b/energy_demand/plotting/_scrap.py index 48747c50..c3586779 100644 --- a/energy_demand/plotting/_scrap.py +++ b/energy_demand/plotting/_scrap.py @@ -1,34 +1,8 @@ -import pandas as pd -import geopandas -from shapely.geometry import Point -import matplotlib.pyplot as plt - -df = pd.DataFrame( - {'City': ['Buenos Aires', 'Brasilia', 'Santiago', 'Bogota', 'Caracas'], - 'Country': ['Argentina', 'Brazil', 'Chile', 'Colombia', 'Venezuela'], - 'Latitude': [-34.58, -15.78, -33.45, 4.60, 10.48], - 'Longitude': [-58.66, -47.91, -70.66, -74.08, -66.86]}) - -df['Coordinates'] = list(zip(df.Longitude, df.Latitude)) - - -df['Coordinates'] = df['Coordinates'].apply(Point) - -gdf = geopandas.GeoDataFrame(df, geometry='Coordinates') -world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres')) - - -# We restrict to South America. -#ax = world[world.continent == 'South America'].plot( -# color='white', edgecolor='black') -ax = world[world.name == "United Kingdom"].plot( - color='white', edgecolor='black') -#uk = world[world.name == "United Kingdom"] - -# We can now plot our GeoDataFrame. -gdf.plot(ax=ax, color='red') - -plt.show() +import math +import numpy as np +import matplotlib.pyplot as plt +import pylab +from scipy.interpolate import interp1d +from scipy.interpolate import spline -print("aa") diff --git a/energy_demand/plotting/basic_plot_functions.py b/energy_demand/plotting/basic_plot_functions.py index bdd49576..828f1abb 100644 --- a/energy_demand/plotting/basic_plot_functions.py +++ b/energy_demand/plotting/basic_plot_functions.py @@ -19,54 +19,93 @@ def cm2inch(*tupl): else: return tuple(i/inch for i in tupl) -def smooth_data(x_list, y_list, num=500, spider=False): +def simple_smooth(x_list, y_list, num=500, spider=False, interpol_kind='quadratic'): + + nr_x_values = len(x_list) + min_x_val = min(x_list) + max_x_val = max(x_list) + + x_values = np.linspace(min_x_val, max_x_val, num=nr_x_values, endpoint=True) + + f2 = interp1d(x_values, y_list, kind=interpol_kind) + + smoothed_data_x = np.linspace( + min_x_val, + max_x_val, + num=num, + endpoint=True) + + smoothed_data_y = f2(smoothed_data_x) + + return smoothed_data_x, smoothed_data_y + +def smooth_data( + x_list, + y_list, + num=500, + spider=False, + interpol_kind='quadratic' + ): """Smooth data x_list : list List with x values y_list : list List with y values + num : int + New number of interpolation points spider : bool Criteria whether spider plot or not - # https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html + + Note: ------ - needs at least 4 entries in lists + - https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html + """ + if spider: - nr_x_values = len(x_list) min_x_val = min(x_list) max_x_val = math.pi * 2 #max is tow pi - x_values = np.linspace(min_x_val, max_x_val, num=nr_x_values, endpoint=True) - - f2 = interp1d(x_values, y_list, kind='quadratic') #quadratic cubic - - smoothed_data_x = np.linspace( + x_values = np.linspace( min_x_val, max_x_val, - num=num, + num=len(x_list), endpoint=True) - else: - nr_x_values = len(x_list) - min_x_val = min(x_list) - max_x_val = max(x_list) - - x_values = np.linspace(min_x_val, max_x_val, num=nr_x_values, endpoint=True) - f2 = interp1d(x_values, y_list, kind='cubic') + f2 = interp1d( + x_values, + y_list, + kind='quadratic') #quadratic cubic - smoothed_data_x = np.linspace( + x_smooth = np.linspace( min_x_val, max_x_val, num=num, endpoint=True) - smoothed_data_y = f2(smoothed_data_x) + #y_smooth = f2(x_smooth) - return smoothed_data_x, smoothed_data_y + else: + # Smooth x data + x_smooth = np.linspace( + min(x_list), + max(x_list), + num=num) + + f2 = interp1d( + x_list, + y_list, + kind=interpol_kind) + + # smooth + y_smooth = f2(x_smooth) + + return x_smooth, y_smooth def smooth_line( input_x_line_data, @@ -78,6 +117,9 @@ def smooth_line( nr_line_points : int represents number of points to make between input_line_data.min and T.max """ + input_x_line_data = np.array(input_x_line_data) + input_y_line_data = np.array(input_y_line_data) + smooth_x_line_data = np.linspace( input_x_line_data.min(), input_x_line_data.max(), diff --git a/energy_demand/plotting/fig_load_profile_dh_multiple.py b/energy_demand/plotting/fig_load_profile_dh_multiple.py index f007fb48..0d233ce3 100644 --- a/energy_demand/plotting/fig_load_profile_dh_multiple.py +++ b/energy_demand/plotting/fig_load_profile_dh_multiple.py @@ -36,10 +36,8 @@ def run( """ fig = plt.figure( figsize=basic_plot_functions.cm2inch(14, 25)) - - ax = fig.add_subplot( - nrows=4, - ncols=2) + + ax = fig.add_subplot(nrows=4, ncols=2) plot_nr = 0 row = -1 @@ -232,8 +230,6 @@ def run( round(std_dev_abs, 2), fontdict=font_additional_info)) plt.title(title_info, loc='left', fontdict=font_additional_info) - #plt.ylabel("hours") - #plt.ylabel("average electricity [GW]") # ------------ # Plot legend @@ -358,14 +354,25 @@ def plot_radar_plot( linestyle='--', linewidth=0.5) + # plot data points + '''ax.scatter(angles, values, color="blue", s=10)''' + + '''line_zeros = np.zeros(len(angles_smoothed)) + ax.fill_between( + angles_smoothed, #zero coordinates + line_zeros, # line 1 + values_smoothed, # line 2 + color='blue', + alpha=0.1)''' + + # Fill below ax.fill( - angles, - values, + angles_smoothed, #angles, + values_smoothed, #values, 'blue', #b alpha=0.1) # Save fig - print("fig_name: " + str(fig_name)) plt.savefig(fig_name) if plotshow: diff --git a/energy_demand/plotting/fig_total_demand_peak.py b/energy_demand/plotting/fig_total_demand_peak.py index 9026a741..7d5f1843 100644 --- a/energy_demand/plotting/fig_total_demand_peak.py +++ b/energy_demand/plotting/fig_total_demand_peak.py @@ -1,6 +1,6 @@ """Plot peak and total demand for """ - +import os import numpy as np import matplotlib.pyplot as plt import pandas as pd @@ -8,6 +8,7 @@ from energy_demand.basic import conversions from energy_demand.technologies import tech_related from energy_demand.plotting import basic_plot_functions +from energy_demand.read_write import write_data def run( data_input, @@ -16,6 +17,8 @@ def run( ): """Plot peak demand and total demand over time in same plot """ + statistics_to_print = [] + # Select period and fueltype fueltype_int = tech_related.get_fueltype_int(fueltype_str) @@ -28,8 +31,9 @@ def run( weather_yrs_peak_demand = [] nr_weather_yrs = list(data_input.keys()) - - print("Weather yrs: " + str(nr_weather_yrs), flush=True) + statistics_to_print.append("_____________________________") + statistics_to_print.append("Weather years") + statistics_to_print.append(str(data_input.keys())) # Iterate weather years for weather_yr, data_weather_yr in data_input.items(): @@ -60,6 +64,10 @@ def run( weather_yrs_total_demand = np.array(weather_yrs_total_demand) weather_yrs_peak_demand = np.array(weather_yrs_peak_demand) + # Calculate std per simulation year + std_total_demand = list(np.std(weather_yrs_total_demand, axis=0)) # across columns calculate std + std_peak_demand = list(np.std(weather_yrs_peak_demand, axis=0)) # across columns calculate std + # Create dataframe if len(nr_weather_yrs) > 2: @@ -114,6 +122,7 @@ def run( period_h_smoothed, df_total_demand_q_05_smoothed = basic_plot_functions.smooth_data(columns, df_total_demand_q_05, num=40000) period_h_smoothed, df_peak_q_95_smoothed = basic_plot_functions.smooth_data(list(columns), df_peak_q_95, num=40000) period_h_smoothed, df_peak_q_05_smoothed = basic_plot_functions.smooth_data(columns, df_peak_q_05, num=40000) + period_h_smoothed, df_peak_2015_smoothed = basic_plot_functions.smooth_data(columns, df_peak_2015, num=40000) except: period_h_smoothed = columns df_total_demand_q_95_smoothed = df_total_demand_q_95 @@ -121,75 +130,143 @@ def run( df_peak_q_95_smoothed = df_peak_q_95 df_peak_q_05_smoothed = df_peak_q_05 tot_demand_twh_2015_smoothed = tot_demand_twh_2015 + df_peak_2015_smoothed = df_peak_2015 else: try: period_h_smoothed, tot_demand_twh_2015_smoothed = basic_plot_functions.smooth_data(columns, tot_demand_twh_2015, num=40000) + period_h_smoothed, df_peak_2015_smoothed = basic_plot_functions.smooth_data(columns, df_peak_2015, num=40000) except: period_h_smoothed = columns tot_demand_twh_2015_smoothed = tot_demand_twh_2015 - + df_peak_2015_smoothed = df_peak_2015 + # -------------- # Two axis figure # -------------- - fig, ax1 = plt.subplots() #figsize=basic_plot_functions.cm2inch(10,10)) + fig, ax1 = plt.subplots( + figsize=basic_plot_functions.cm2inch(15, 10)) + ax2 = ax1.twinx() # Axis label - ax1.set_xlabel('years') - ax2.set_ylabel('Peak hour demand (GW)', color='blue') - ax1.set_ylabel('Total demand (TWh)', color='tomato') - - # Make the y-axis label, ticks and tick labels match the line color. - ax1.tick_params('y', colors='tomato') - ax2.tick_params('y', colors='blue') - + ax1.set_xlabel('Years') + ax2.set_ylabel('Peak hour {} demand (GW)'.format(fueltype_str), color='black') + ax1.set_ylabel('Total {} demand (TWh)'.format(fueltype_str), color='black') + + # Make the y-axis label, ticks and tick labels match the line color.ยจ + color_axis1 = 'lightgrey' + color_axis2 = 'blue' + + ax1.tick_params('y', colors='black') + ax2.tick_params('y', colors='black') + if len(nr_weather_yrs) > 2: - + # ----------------- # Uncertainty range total demand # ----------------- - ax1.plot(period_h_smoothed, df_total_demand_q_05_smoothed, color='tomato', linestyle='--', linewidth=0.5, label="0.05_total_demand") - ax1.plot(period_h_smoothed, df_total_demand_q_95_smoothed, color='tomato', linestyle='--', linewidth=0.5, label="0.95_total_demand") + '''ax1.plot( + period_h_smoothed, + df_total_demand_q_05_smoothed, + color='tomato', linestyle='--', linewidth=0.5, label="0.05_total_demand")''' + + '''ax1.plot( + period_h_smoothed, + df_total_demand_q_95_smoothed, + color=color_axis1, linestyle='--', linewidth=0.5, label="0.95_total_demand") + ax1.fill_between( period_h_smoothed, #x df_total_demand_q_95_smoothed, #y1 df_total_demand_q_05_smoothed, #y2 alpha=.25, - facecolor="tomato", - label="uncertainty band demand") - + facecolor=color_axis1, + label="uncertainty band demand")''' + # ----------------- # Uncertainty range peaks # ----------------- - ax2.plot(period_h_smoothed, df_peak_q_05_smoothed, color='blue', linestyle='--', linewidth=0.5, label="0.05_peak") - ax2.plot(period_h_smoothed, df_peak_q_95_smoothed, color='blue', linestyle='--', linewidth=0.5, label="0.95_peak") - ax2.fill_between( + ##ax2.plot(period_h_smoothed, df_peak_q_05_smoothed, color=color_axis2, linestyle='--', linewidth=0.5, label="0.05_peak") + ##ax2.plot(period_h_smoothed, df_peak_q_95_smoothed, color=color_axis2, linestyle='--', linewidth=0.5, label="0.95_peak") + ax2.plot( + period_h_smoothed, + df_peak_2015_smoothed, + color=color_axis2, linestyle="--", linewidth=0.4) + + # Error bar of bar charts + ax2.errorbar(columns, df_peak_2015, linewidth=0.5, color='black', yerr=std_peak_demand, linestyle="None") + + # Error bar bar plots + ax1.errorbar( + columns, tot_demand_twh_2015, linewidth=0.5, color='black', yerr=std_total_demand, linestyle="None") + + '''ax2.fill_between( period_h_smoothed, #x df_peak_q_95_smoothed, #y1 df_peak_q_05_smoothed, #y2 alpha=.25, facecolor="blue", - label="uncertainty band peak") + label="uncertainty band peak")''' + + # Total demand bar plots + ##ax1.plot(period_h_smoothed, tot_demand_twh_2015_smoothed, color='tomato', linestyle='-', linewidth=2, label="tot_demand_weather_yr_2015") + ax1.bar( + columns, + tot_demand_twh_2015, + width=2, + alpha=1, + align='center', + color=color_axis1, + label="total {} demand".format(fueltype_str)) + + statistics_to_print.append("_____________________________") + statistics_to_print.append("total demand per model year") + statistics_to_print.append(str(tot_demand_twh_2015)) - # Plot weather year total demand - ax1.plot(period_h_smoothed, tot_demand_twh_2015_smoothed, color='tomato', linestyle='-', linewidth=2, label="tot_demand_weather_yr_2015") + # Line of peak demand + #ax2.plot(columns, df_peak, color=color_axis2, linestyle='--', linewidth=0.5, label="peak_0.95") + ax2.plot( + period_h_smoothed, + df_peak_2015_smoothed, + color=color_axis2, linestyle='-', linewidth=2, label="{} peak demand (base weather yr)".format(fueltype_str)) - # plot peak (all values) - #ax2.plot(columns, df_peak, color='blue', linestyle='--', linewidth=0.5, label="peak_0.95") - ax2.plot(columns, df_peak_2015, color='blue', linestyle='-', linewidth=2, label="peak_weather_yr_2015") + statistics_to_print.append("_____________________________") + statistics_to_print.append("peak demand per model year") + statistics_to_print.append(str(df_peak_2015)) + + # Scatter plots of peak demand + ax2.scatter( + columns, + df_peak_2015, + marker='o', s=20, color=color_axis2, alpha=1) ax1.legend( prop={ 'family':'arial', 'size': 10}, - loc='right', + loc='upper center', + bbox_to_anchor=(0.9, -0.1), frameon=False, shadow=True) + ax2.legend( prop={ 'family':'arial', 'size': 10}, - loc='left', + loc='upper center', + bbox_to_anchor=(0.1, -0.1), frameon=False, shadow=True) - plt.show() + + # More space at bottom + #fig.subplots_adjust(bottom=0.4) + fig.tight_layout() + + plt.savefig(fig_name) + plt.close() + + # Write info to txt + write_data.write_list_to_txt( + os.path.join(fig_name.replace(".pdf", ".txt")), + statistics_to_print) + print("--") diff --git a/energy_demand/plotting/result_mapping.py b/energy_demand/plotting/result_mapping.py index 95bcdd9b..97b96950 100644 --- a/energy_demand/plotting/result_mapping.py +++ b/energy_demand/plotting/result_mapping.py @@ -608,6 +608,7 @@ def plot_lad_national( # Show figure if plotshow: plt.show() + plt.close() def merge_data_to_shp(shp_gdp, merge_data, unique_merge_id): diff --git a/energy_demand/result_processing/scenario_charts.py b/energy_demand/result_processing/scenario_charts.py index 156695e2..976053df 100644 --- a/energy_demand/result_processing/scenario_charts.py +++ b/energy_demand/result_processing/scenario_charts.py @@ -56,8 +56,6 @@ def process_scenarios(path_to_scenarios, year_to_model=2015): # ------------------------------- scenario_data = {} for scenario in scenarios: - - # Add scenario name to folder scenario_data[scenario] = {} path_to_result_files = os.path.join( diff --git a/energy_demand/result_processing/single_scenario.py b/energy_demand/result_processing/single_scenario.py index 229030d7..a05c108c 100644 --- a/energy_demand/result_processing/single_scenario.py +++ b/energy_demand/result_processing/single_scenario.py @@ -6,13 +6,13 @@ from energy_demand.read_write import data_loader, read_data from energy_demand.basic import date_prop from energy_demand.plotting import plotting_results, result_mapping -from energy_demand.basic import logger_setup, basic_functions +from energy_demand.basic import basic_functions from energy_demand.basic import lookup_tables from energy_demand.plotting import fig_weather_variability_priod from energy_demand.plotting import fig_total_demand_peak def main( - path_data_energy_demand, + path_data_ed, path_shapefile_input, plot_crit_dict, base_yr, @@ -22,7 +22,7 @@ def main( Arguments ---------- - path_data_energy_demand : str + path_data_ed : str Path to results path_shapefile_input : str Path to shapefile @@ -34,9 +34,12 @@ def main( Year to generate comparison plots """ print("...Start creating plots") + data = {} + # ------------------------------------------------------------ # Get all yealy results - all_result_folders = os.listdir(path_data_energy_demand) + # ------------------------------------------------------------ + all_result_folders = os.listdir(path_data_ed) weather_yrs = [] for result_folder in all_result_folders: @@ -46,22 +49,21 @@ def main( except ValueError: pass - # ------------------------------ + # ------------------------------------------------------------ # Plotting weather variability results - # ------------------------------ + # ------------------------------------------------------------ if plot_crit_dict['plot_weather_day_year']: # Container to store all data of weather years weather_yr_container = defaultdict(dict) - data = {} data['lookups'] = lookup_tables.basic_lookups() - path_out_plots = os.path.join(path_data_energy_demand, "PDF_weather_varability") + path_out_plots = os.path.join(path_data_ed, "PDF_weather_varability") basic_functions.del_previous_setup(path_out_plots) basic_functions.create_folder(path_out_plots) data['enduses'], data['assumptions'], data['reg_nrs'], data['regions'] = data_loader.load_ini_param( - os.path.join(path_data_energy_demand)) + os.path.join(path_data_ed)) # Other information is read in data['assumptions']['seasons'] = date_prop.get_season(year_to_model=2015) @@ -73,7 +75,7 @@ def main( for weather_yr in weather_yrs: results_container = read_data.read_in_results( - os.path.join(path_data_energy_demand, str(weather_yr), 'model_run_results_txt'), + os.path.join(path_data_ed, str(weather_yr), 'model_run_results_txt'), data['assumptions']['seasons'], data['assumptions']['model_yeardays_daytype']) @@ -84,14 +86,15 @@ def main( weather_yr_container['tot_fueltype_h'][weather_yr] = tot_fueltype_h # -------------------------------------------- - # Plots + # Plot peak demand and total demand per fueltype # -------------------------------------------- - # plot over time total demand and peak - fig_total_demand_peak.run( - data_input=weather_yr_container['tot_fueltype_h'], - fueltype_str='electricity', - fig_name=os.path.join( - path_out_plots, "tot_fueltype_h.pdf")) + for fueltype_str in data['lookups']['fueltypes'].keys(): + # plot over time total demand and peak + fig_total_demand_peak.run( + data_input=weather_yr_container['tot_fueltype_h'], + fueltype_str=fueltype_str, + fig_name=os.path.join( + path_out_plots, "tot_{}_h.pdf".format(fueltype_str))) # plot over period of time across all weather scenario @@ -102,36 +105,28 @@ def main( period_h=list(range(200,500)), #period to plot fig_name=os.path.join( path_out_plots, "weather_var_period.pdf")) - - - - - - else: pass - # Execute script to generate PDF results + # ------------------------------------------------------------ + # Calculate results for every weather year + # ------------------------------------------------------------ for weather_yr in weather_yrs: - path_data_energy_demand_weather_yr = os.path.join(path_data_energy_demand, str(weather_yr)) + path_data_weather_yr = os.path.join(path_data_ed, str(weather_yr)) - # Set up logger - #logger_setup.set_up_logger( - # os.path.join( - # path_data_energy_demand_weather_yr, "plotting.log")) # Simulation information is read in from .ini file for results data['enduses'], data['assumptions'], data['reg_nrs'], data['regions'] = data_loader.load_ini_param( - os.path.join(path_data_energy_demand)) + os.path.join(path_data_ed)) # ------------------ # Load necessary inputs for read in # ------------------ data = {} data['local_paths'] = data_loader.get_local_paths( - path_data_energy_demand_weather_yr) + path_data_weather_yr) data['result_paths'] = data_loader.get_result_paths( - os.path.join(path_data_energy_demand_weather_yr)) + os.path.join(path_data_weather_yr)) data['lookups'] = lookup_tables.basic_lookups() # --------------- @@ -145,17 +140,15 @@ def main( # Simulation information is read in from .ini file for results data['enduses'], data['assumptions'], data['reg_nrs'], data['regions'] = data_loader.load_ini_param( - os.path.join(path_data_energy_demand)) + os.path.join(path_data_ed)) # Other information is read in data['assumptions']['seasons'] = date_prop.get_season(year_to_model=2015) data['assumptions']['model_yeardays_daytype'], data['assumptions']['yeardays_month'], data['assumptions']['yeardays_month_days'] = date_prop.get_yeardays_daytype(year_to_model=2015) - # Read scenario data data['scenario_data'] = {} - data['scenario_data']['population'] = read_data.read_scenaric_population_data( - os.path.join(path_data_energy_demand, 'model_run_pop')) + os.path.join(path_data_ed, 'model_run_pop')) # -------------------------------------------- # Reading in results from different model runs