- added that reg_names are passed along to ini file

- write function to append  model results to shapefile

energy_demand/basic/basic_functions.py

@@ -17,6 +17,28 @@ def rmse(predictions, targets):
     """
     return np.sqrt(((predictions - targets) ** 2).mean())
 
+def array_to_dict(result_array, lu_reg):
+    """Convert an array with regions to dict
+    with region as key
+
+    Arguments
+    ---------
+    result_array : array
+        Results in region_array
+    lu_reg : list
+        List with all regions (oder is the same)
+
+    Returns
+    --------
+    result_dict : dict
+        reg, value
+    """
+    result_dict = {}
+    for reg_array_nr, region_name in enumerate(lu_reg):
+        result_dict[region_name] = result_array[reg_array_nr]
+
+    return result_dict
+
 def create_folder(path_folder, name_subfolder=None):
     """Creates folder or subfolder
 
@@ -71,7 +93,7 @@ def del_previous_setup(path_folder):
     if os.path.exists(path_folder):
         shutil.rmtree(path_folder, ignore_errors=True)
         logging.info("... deleted previous scenario results")
-        print("...Deleted previous scenario results")
+        print("... deleted previous scenario results")
     else:
         pass
 
@@ -153,4 +175,4 @@ def get_month_from_int(month_int):
     elif month_int == 12:
         month_str = 'Dec'
 
-    return str(month_str)
+    return str(month_str)

energy_demand/charts/figure_HHD_gas_demand.py

@@ -8,12 +8,10 @@
 Note: This could be improved by calculating the HDD for every region and multiplying
 by the population. This would take into consideration the different population
 """
-import os
 import logging
 import numpy as np
 from scipy import stats
 import matplotlib.pyplot as plt
-from energy_demand.profiles import hdd_cdd
 from energy_demand.scripts import s_raw_weather_data
 from energy_demand.plotting import plotting_program
 from energy_demand.geography.weather_station_location import get_closest_station
@@ -517,8 +515,8 @@ def lin_func(x, slope, intercept):
     font_additional_info = {'family': 'arial', 'color': 'black', 'weight': 'normal', 'size': 8}
     plt.xlabel("National daily gas demand [GWh per day]")
     plt.ylabel("Heating degree days [mio]")
-    plt.title("Correlation between national gas demand and hdd (r_value:  {}".format(r_value), fontdict=font_additional_info)
-    plt.text(1, 0.55, "y = {} x + {}".format(slope, intercept), fontdict=font_additional_info)
+    plt.title("LAD HDD * LADPOP vs Gas_y_uk (r_value: {}".format(round(r_value, 3)), fontdict=font_additional_info)
+    #plt.text(10, 10, "y = {} x + {}".format(round(slope, 2), round(intercept, 2)), fontdict=font_additional_info)
 
     # Tight layout
     #plt.tight_layout()

energy_demand/energy_model.py

@@ -16,6 +16,7 @@
 from energy_demand.profiles import load_profile
 from energy_demand.initalisations import helpers
 from energy_demand.profiles import load_factors as lf
+from energy_demand.charts import figure_HHD_gas_demand
 
 class EnergyModel(object):
     """EnergyModel of a simulation yearly run
@@ -196,18 +197,17 @@ def __init__(self, region_names, data):
         self.reg_load_factor_yd = reg_load_factor_yd
         self.reg_load_factor_seasons = reg_load_factor_seasons
 
-        #-------------------
+        # ------------------------------
         # TESTING
-        #-------------------
+        # ------------------------------
         testing.test_region_selection(self.ed_fueltype_regs_yh)
 
-        # ---
-        # Other charts
-        # ---
-        plot_HDD_chart = True
+        # ------------------------------
+        # Chart HDD * Pop vs actual gas demand
+        # ------------------------------
+        plot_HDD_chart = False
         if plot_HDD_chart == True:
             logging.info("plot figure HDD comparison")
-            from energy_demand.charts import figure_HHD_gas_demand
             figure_HHD_gas_demand.main(region_names, weather_regions, data)
 
 def simulate_region(region_name, data, weather_regions):

energy_demand/geography/write_shp.py

@@ -0,0 +1,89 @@
+"""Append result to shapefile
+"""
+import logging
+import shapefile
+
+def write_result_shapefile(lad_geometry_shp, out_shape, field_names, csv_results):
+    """
+    Join result attributes to LAD geography with
+    pyhape library
+
+    Arguments
+    ---------
+    lad_geometry_shp : str
+        Path to LAD shapefile
+    out_shape : str
+        Path to new shapefile
+    field_names : list
+        list with new attribute field name
+    csv_results : list
+        list with result dicts
+
+    Info
+    -----
+    pip install pyshp
+    https://github.com/GeospatialPython/pyshp#reading-shapefiles-from-file-like-objects
+
+    http://www.qgistutorials.com/en/docs/performing_table_joins_pyqgis.html
+    """
+    # Read in our existing shapefile
+    lad_geometry_shp_name = lad_geometry_shp[:-3]
+    myshp = open(lad_geometry_shp_name + "shp", "rb")
+    mydbf = open(lad_geometry_shp_name + "dbf", "rb")
+
+    r = shapefile.Reader(shp=myshp, dbf=mydbf)
+
+    # Create a new shapefile in memory
+    w = shapefile.Writer()
+
+    # Copy over the existing fields
+    w.fields = list(r.fields)
+
+    # --------------
+    # Add new fields
+    # --------------
+    for field_name in field_names:
+        w.field(field_name, "F", decimal=10) #Float
+
+    # Get position of field 'geo_code'
+    position = 0
+    for field_name in r.fields[1:]:
+        if field_name[0] == 'geo_code':
+            position_field_name = position
+        else:
+            position += 1
+
+    # --------------------------
+    # Join fields programatically
+    # --------------------------
+    # Loop through each record, add a column and get results
+    for rec in r.records():
+
+        # Get geocode for row
+        geo_code = rec[position_field_name]
+
+        # Iterate result entries in list
+        for result_per_field in csv_results:
+
+            # Iterate result entries and add
+            #for _ in result_per_field:
+            try:
+                result_csv = result_per_field[geo_code]
+            except KeyError:
+                # No results
+                result_csv = 0
+                logging.warning(
+                    "No value found for region '%s' in merging shapefile results", geo_code)
+
+                rec.append(result_csv)
+
+        # Add the modified record to the new shapefile
+        w.records.append(rec)
+
+    # Copy over the geometry without any changes
+    w._shapes.extend(r.shapes())
+
+    # Save as a new shapefile (or write over the old one)
+    w.save(out_shape)
+    logging.info("... finished writing shp")
+    return

energy_demand/main.py

@@ -118,15 +118,13 @@ def energy_demand_model(data, fuel_in=0, fuel_in_elec=0):
     # Parameters not defined within smif
     data['assumptions'] = non_param_assumptions.load_non_param_assump(
         data['sim_param']['base_yr'], data['paths'], data['enduses'], data['lookups'], data['fuels'])
-    
+
     # Parameters defined within smif
     param_assumptions.load_param_assump(data['paths'], data['assumptions'])
-    #data['assumptions'] = read_data.read_param_yaml(data['paths']['yaml_parameters'])
 
     data['assumptions']['seasons'] = date_prop.read_season(year_to_model=2015)
     data['assumptions']['model_yeardays_daytype'], data['assumptions']['yeardays_month'], data['assumptions']['yeardays_month_days'] = date_prop.get_model_yeardays_datype(year_to_model=2015)
 
-
     data['tech_lp'] = data_loader.load_data_profiles(data['paths'], data['local_paths'], data['assumptions']['model_yeardays'], data['assumptions']['model_yeardays_daytype'])
     data['assumptions']['technologies'] = non_param_assumptions.update_assumptions(data['assumptions']['technologies'], data['assumptions']['eff_achiev_f']['factor_achieved'])
     data['weather_stations'], data['temp_data'] = data_loader.load_temp_data(data['local_paths'])
@@ -172,7 +170,8 @@ def energy_demand_model(data, fuel_in=0, fuel_in_elec=0):
         data['sim_param'],
         data['enduses'],
         data['assumptions'],
-        data['reg_nrs'])
+        data['reg_nrs'],
+        data['lu_reg'])
 
     for sim_yr in data['sim_param']['simulated_yrs']:
         data['sim_param']['curr_yr'] = sim_yr
@@ -247,7 +246,6 @@ def energy_demand_model(data, fuel_in=0, fuel_in_elec=0):
                 model_run_object.tot_peak_enduses_fueltype)
 
 
-
     # ------------------
     # Load necessary inputs for read in
     # ------------------
@@ -257,13 +255,12 @@ def energy_demand_model(data, fuel_in=0, fuel_in_elec=0):
     data['lookups'] = data_loader.load_basic_lookups()
 
     # Simulation information is read in from .ini file for results
-    data['sim_param'], data['enduses'], data['assumptions'], data['reg_nrs'] = data_loader.load_sim_param_ini(data['local_paths']['data_results'])
+    data['sim_param'], data['enduses'], data['assumptions'], data['reg_nrs'], data['lu_reg'] = data_loader.load_sim_param_ini(data['local_paths']['data_results'])
 
     # Other information is read in
     data['assumptions']['seasons'] = date_prop.read_season(year_to_model=2015)
     data['assumptions']['model_yeardays_daytype'], data['assumptions']['yeardays_month'], data['assumptions']['yeardays_month_days'] = date_prop.get_model_yeardays_datype(year_to_model=2015)
 
-
     # --------------------------------------------
     # Reading in results from different model runs
     # --------------------------------------------
@@ -285,4 +282,4 @@ def energy_demand_model(data, fuel_in=0, fuel_in_elec=0):
         data['sim_param'],
         data['enduses'])
 
-    logging.info("... Finished running Energy Demand Model")
+    logging.info("... Finished running Energy Demand Model")

energy_demand/plotting/plotting_results.py

@@ -833,7 +833,6 @@ def plt_fuels_peak_h(tot_fuel_dh_peak, lookups, path_plot_fig):
     # -
     # Axis
     # -
-
     base_yr = 2015
     major_interval = 10
     minor_interval = 5

energy_demand/read_and_plot_results.py

@@ -1,9 +1,11 @@
 """Read in model results and plot results
 """
+import os
 import logging
-from energy_demand.read_write import data_loader, read_data
+from energy_demand.read_write import data_loader, read_data, write_data
 from energy_demand.basic import date_prop
 from energy_demand.plotting import plotting_results
+from energy_demand.basic import logger_setup
 
 def main(local_data_path):
     """Read in all results and plot PDFs
@@ -13,6 +15,7 @@ def main(local_data_path):
     local_data_path : str
         Path to results
     """
+    logger_setup.set_up_logger(os.path.join(local_data_path, "logging_plotting.log"))
     # ------------------
     # Load necessary inputs for read in
     # ------------------
@@ -21,7 +24,7 @@ def main(local_data_path):
     data['lookups'] = data_loader.load_basic_lookups()
 
     # Simulation information is read in from .ini file for results
-    data['sim_param'], data['enduses'], data['assumptions'], data['reg_nrs'] = data_loader.load_sim_param_ini(
+    data['sim_param'], data['enduses'], data['assumptions'], data['reg_nrs'], data['lu_reg'] = data_loader.load_sim_param_ini(
         data['local_paths']['data_results'])
 
     # Other information is read in
@@ -37,6 +40,15 @@ def main(local_data_path):
         data['assumptions']['seasons'],
         data['assumptions']['model_yeardays_daytype'])
 
+    # ----------------
+    # Write results to CSV files and merge with shapefile
+    # ----------------
+    write_data.create_shp_results(
+        results_container,
+        data['local_paths'],
+        data['lookups'],
+        data['lu_reg'])
+
     # ------------------------------
     # Plotting results
     # ------------------------------

energy_demand/read_write/data_loader.py

@@ -20,7 +20,7 @@ def load_sim_param_ini(path):
     config.read(os.path.join(path, 'model_run_sim_param.ini'))
 
     reg_nrs = int(config['SIM_PARAM']['reg_nrs'])
-
+    lu_reg = ast.literal_eval(config['REGIONS']['lu_reg'])
     sim_param = {}
     sim_param['base_yr'] = int(config['SIM_PARAM']['base_yr'])
     sim_param['simulated_yrs'] = ast.literal_eval(config['SIM_PARAM']['simulated_yrs'])
@@ -37,7 +37,7 @@ def load_sim_param_ini(path):
     enduses['ss_all_enduses'] = ast.literal_eval(config['ENDUSES']['ss_all_enduses'])
     enduses['is_all_enduses'] = ast.literal_eval(config['ENDUSES']['is_all_enduses'])
 
-    return sim_param, enduses, assumptions, reg_nrs
+    return sim_param, enduses, assumptions, reg_nrs, lu_reg
 
 def read_national_real_elec_data(path_to_csv):
     """Read in national consumption from csv file
@@ -277,14 +277,16 @@ def load_local_paths(path):
             path, '_processed_data'),
         'data_results': os.path.join(
             path, '_result_data'),
+        'lad_shapefile_2011': os.path.join(
+            path, '_raw_data', 'C_LAD_geography', 'infuse_dist_lyr_2011.shp'),
         'path_post_installation_data': os.path.join(
             path, '_processed_data', '_post_installation_data'),
         'data_processed_disaggregated': os.path.join(
             path, '_processed_data', '_post_installation_data', 'disaggregated'),
         'path_sigmoid_data': os.path.join(
             path, '_processed_data', 'sigmoid_data'),
         'data_results_model_runs': os.path.join(
-            path, '_result_data', "model_run_results_txt"),
+            path, '_result_data', 'model_run_results_txt'),
         'dir_changed_weather_station_data': os.path.join(
             path, '_processed_data', '_post_installation_data', 'weather_station_data'),
         'changed_weather_station_data': os.path.join(
@@ -306,7 +308,9 @@ def load_local_paths(path):
         'dir_services': os.path.join(
             path, '_processed_data', 'services'),
         'data_results_PDF': os.path.join(
-            path, '_result_data', 'PDF')}
+            path, '_result_data', 'PDF'),
+        'data_results_shapefiles': os.path.join(
+            path, '_result_data', 'result_shapefiles')}
 
     return paths
 

energy_demand/read_write/read_data.py

@@ -34,6 +34,7 @@ def read_in_results(path_runs, lookups, seasons, model_yeardays_daytype):
     # -------------
     # Fuels
     # -------------
+    logging.info("... Reading in fuels")
     results_container['results_every_year'] = read_results_yh(
         lookups['fueltypes_nr'], path_runs)
 
@@ -46,6 +47,7 @@ def read_in_results(path_runs, lookups, seasons, model_yeardays_daytype):
     # -------------
     # Load factors
     # -------------
+    logging.info("... Reading in load factors")
     results_container['load_factors_y'] = read_lf_y(
         os.path.join(path_runs, "result_reg_load_factor_y"))
     results_container['load_factors_yh'] = read_lf_y(
@@ -64,6 +66,7 @@ def read_in_results(path_runs, lookups, seasons, model_yeardays_daytype):
     # -------------
     # Post-calculations
     # -------------
+    logging.info("... generating post calculations with read results")
     # Calculate average per season and fueltype for every fueltype
     av_season_daytype_cy = {}
     season_daytype_cy = {}
@@ -130,6 +133,7 @@ def read_results_yh(fueltypes_nr, path_to_folder):
 
     # Iterate files in folder
     for file_path in all_txt_files_in_folder:
+        logging.info("... file_path: " + str(file_path))
         try:
             path_file_to_read = os.path.join(path_to_folder, file_path)
             file_path_split = file_path.split("__")