<jupyter_start><jupyter_text>CO2 calculation bottom up method - From power plant data to country aggregated CO2 intensity of electricity generation - In this script we calculated an CO2 emission factor per country out of hourly power plant generation and yearly published emission data on plant level.The used method follows the idea to calculated a CI for specific power plants. In a second step an representative sample of power plants for a country is build and an CI for each technology and country is calculated.  Script setup<jupyter_code>import os
import logging

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
from IPython.display import Image 

%matplotlib inline
plt.style.use('seaborn')
plt.rcParams['figure.figsize'] = [15, 10]

#helpers
def change_ProductionTypeName (ProductionTypeName):
    return ProductionTypeName.replace({
            'Fossil Hard coal': 'hard_coal', 
            'Fossil Brown coal/Lignite': 'lignite', 
            'Fossil Gas': 'gas',
            'Fossil Coal-derived gas': 'other_fossil',
            'Fossil Peat': 'other_fossil',
            'Nuclear': 'nuclear',
            'Biomass': 'biomass',
            'Waste': 'waste',
            'Geothermal': 'geothermal',
            'Marine': 'marine',
            'Other': 'other_fossil',
            'Hydro Pumped Storage': 'hydro',
            'Hydro Run-of-river and poundage': 'hydro',
            'Hydro Water Reservoir': 'hydro',
            'Fossil Oil': 'oil',
            'Fossil Oil shale': 'oil', 
            'Solar': 'solar',
            'Wind Onshore': 'wind_onshore',
            'Wind Offshore': 'wind_offshore',
            'Other renewable': 'other_renewable'}, regex = True, inplace = False)<jupyter_output><empty_output><jupyter_text>Data directory preparention Create input, processed and output folders if they don't existIf the paths are relative, the corresponding folders will be created inside the current working directory.<jupyter_code>input_directory_path = os.path.join('input')
Bootom_up_methode_input_directory_path = os.path.join('input', 'Bootom_up_methode')
processed_directory_path = 'processed'
output_directory_path = os.path.join('output')

os.makedirs(input_directory_path, exist_ok=True)
os.makedirs(Bootom_up_methode_input_directory_path, exist_ok=True)
os.makedirs(processed_directory_path, exist_ok=True)
os.makedirs(output_directory_path, exist_ok=True)<jupyter_output><empty_output><jupyter_text>Data file preperation The directory `input/Bootom_up_method` should contain all necessary raw data files.Based on the developed method the following data sets needed:1) Matcher for power plants in Entso e and EUTL Data for EuropeThe matching was performed by a manual process (see description).File - > Matching_Entso_EUTL_EU.csv2) EU Emissions Data (EUTL)Data are provided in the report section. See the following link:https://ec.europa.eu/clima/policies/ets/registry_entab-0-1For example: ->Documentation->Reports->Verified Emissions for 2019File - > "Verified Emissions for 2019" -> verified_emissions_2018_en.xlsx -> converted to .csvAnother way to check the data is directly through the European Union Transaction Log:https://ec.europa.eu/clima/ets/napMgt.do?languageCode=en3) ENTSO-E production per unit Production per Unit from ENTSO-EOriginal data: ENTSO-E Transparency Platform, Actual Generation per Generation Unit Available online: https://transparency.entsoe.eu/generation/r2/actualGenerationPerGenerationUnit/show (accessed on Apr 29, 2020).Processed with the following script that converts the data to hourly resolution (private script from INATECH):https://github.com/INATECH-CIG/entso-e_GenerationOutputPerUnitFile - > gen_data.csv -> hourly generation data per unit\File - > unit_data.csv -> information about the generation units4) ENTSO-E production per typeProduction per Unit from ENTSO-EOriginal data: ENTSO-E Transparency Platform, Actual Generation per Production Type Available online: https://transparency.entsoe.eu/generation/r2/actualGenerationPerProductionType/showProcessed with the following script that converts the data to hourly resolution (private script from INATECH):.... need to be updated .... File - ><jupyter_code># Checks if the the input directories are empty or not
# Checks all filenames in the input directory

if not os.listdir(Bootom_up_methode_input_directory_path) :
    print("The directory for the bootom up method is empty. Please provide the data to the directory as described in the instructions above.")


filenames = [os.path.join(Bootom_up_methode_input_directory_path, fn) for fn in os.listdir(Bootom_up_methode_input_directory_path)]

print(filenames)<jupyter_output>['input\\Bootom_up_methode\\entsoe_unit_gen_data.csv', 'input\\Bootom_up_methode\\gen_data.csv', 'input\\Bootom_up_methode\\matching_ENTSOE_EU_ETS.csv', 'input\\Bootom_up_methode\\Matching_Entso_EUTL_EU.csv', 'input\\Bootom_up_methode\\Matching_Entso_EUTL_EU_2.csv', 'input\\Bootom_up_methode\\Matching_Entso_EUTL_LCPD.csv', 'input\\Bootom_up_methode\\nrg_cb_pem.csv', 'input\\Bootom_up_methode\\time_series_60min_multiindex.csv', 'input\\Bootom_up_methode\\unit_data.csv', 'input\\Bootom_up_methode\\verified_emissions_2018_en.csv', 'input\\Bootom_up_methode\\verified_emissions_2018_en.xlsx']<jupyter_text>Load data functions<jupyter_code>def load_matching_data_EU(path, fn):
    """
    Matching List for EU power plants with Entso e identifier and the EUTL identifier.
        
    Parameters
    ----------
    path: str
        path to data
    fn : str
        filename
        
    """
    
    logging.info(f'Loading data from {fn}')
    
    df = pd.read_csv(os.path.join(path, fn), sep = ',', header = 0, index_col=0)

    return df

def load_EUTL_data(path, fn):
    """
    EU Emissions Data (EUTL)
        
    Parameters
    ----------
    path: str
        path to data
    fn : str
        filename
        
    """
    
    logging.info(f'Loading data from {fn}')
    
    df = pd.read_csv(os.path.join(path, fn),sep = ';', header = 13, encoding='ISO-8859-1')

    return df

def load_generation_data(path, fn):
    """
    Entso e gernation data per unit
        
    Parameters
    ----------
    path: str
        path to data
    fn : str
        filename
        
    """
    
    logging.info(f'Loading data from {fn}')
    
    df = pd.read_csv(os.path.join(path, fn),sep = ',',index_col=0,parse_dates=True)
    

    return df

def load_unit_data(path, fn):
    """
    Entso e generation unit information
        
    Parameters
    ----------
    path: str
        path to data
    fn : str
        filename
        
    """
    
    logging.info(f'Loading data from {fn}')
    
    df = pd.read_csv(os.path.join(path, fn),sep = ',',index_col=0)
    
    # Rename production type name according to own convention
    df.ProductionTypeName = change_ProductionTypeName(df.ProductionTypeName)
    
    # set name for the index
    df.index.set_names('GenerationUnitEIC', inplace=True)

    return df

def load_timeseries_opsd(path, fn):
    """
    Read data from OPSD time-series package own modification.

    
    """
    generation = (pd.read_csv(os.path.join(path, fn), index_col=[0], header=[0, 1, 2, 3, 4, 5], parse_dates=True)
                    .dropna(how="all", axis=0))
        
  
    
    generation = generation.rename(columns={'GB_UKM' : 'GB'})
    generation = generation.drop(columns='cet_cest_timestamp', level='region')
       
    
    return generation<jupyter_output><empty_output><jupyter_text>Load data sets  Load power plant generation data.<jupyter_code>generation_per_unit = load_generation_data(Bootom_up_methode_input_directory_path, 'gen_data.csv')<jupyter_output><empty_output><jupyter_text>Load power plant unit inforamtion (capacity, name, etc.)<jupyter_code>generation_unit_info = load_unit_data(Bootom_up_methode_input_directory_path, 'unit_data.csv')<jupyter_output><empty_output><jupyter_text>Load CO2_emissions data from EUTL dataset<jupyter_code>EUTL_emissions = load_EUTL_data(Bootom_up_methode_input_directory_path, 'verified_emissions_2018_en.csv')<jupyter_output><empty_output><jupyter_text>Load machting information for power plant data<jupyter_code>unit_matching_EU = load_matching_data_EU(Bootom_up_methode_input_directory_path, 'matching_ENTSOE_EU_ETS.csv')<jupyter_output><empty_output><jupyter_text>Method  Yearly power generation per unit Store the yearly generation per power plant to the power plant information data set.<jupyter_code>generation_unit_info['generation_2018'] = generation_per_unit.sum()
generation_unit_info[generation_unit_info.PowerSystemResourceName == 'Theiß 2']<jupyter_output><empty_output><jupyter_text>Matching ENTSO-E ID and EUTL ID Combining matching information (ENTSO-E ID and EUTL ID) to the generation unit info. Keep only matched information.<jupyter_code>generation_unit_info_matched = pd.merge(generation_unit_info, unit_matching_EU, left_on='GenerationUnitEIC', right_on='eic_g', how='inner')<jupyter_output><empty_output><jupyter_text>Function for connecting emissions and production data<jupyter_code>def connect_CO2(df, ETS):
    # caring for characteristics of the dataset: EUTL-ID IS NOT UNIQUE - only countrywise.
    # removing power plants from the dataset which are not covered in ETS:
    miss_country = []
    for i in set(df.EUTL_countrycode):
        if sum(i == ETS['REGISTRY_CODE']) == 0:
            print(i)
            miss_country.append(i)
            #removing this power plants from the list:
    df = df.query('EUTL_countrycode not in @miss_country')
    
    # removing power plants where match could not be found:
    miss_match = []
    for j in df.eic_g:
        row = df[df.eic_g == j]
        if len(ETS.query('REGISTRY_CODE == @row.EUTL_countrycode.iloc[0]')\
                  .query('INSTALLATION_IDENTIFIER == @row.EUTL_ID.iloc[0]')) == 0:
            print(j)
            miss_match.append(j)
            
    df = df.query('eic_g not in @miss_match')
    # apply matching:
    df['verified_emissions_18'] = df.apply(lambda x: ETS.query('REGISTRY_CODE==@x.EUTL_countrycode')\
                                    .query('INSTALLATION_IDENTIFIER==@x.EUTL_ID')['VERIFIED_EMISSIONS_2018']\
                                    .iloc[0], axis = 1)
    df['verified_emissions_17'] = df.apply(lambda x: ETS.query('REGISTRY_CODE==@x.EUTL_countrycode')\
                                    .query('INSTALLATION_IDENTIFIER==@x.EUTL_ID')['VERIFIED_EMISSIONS_2017']\
                                    .iloc[0], axis = 1)
    df['verified_emissions_16'] = df.apply(lambda x: ETS.query('REGISTRY_CODE==@x.EUTL_countrycode')\
                                    .query('INSTALLATION_IDENTIFIER==@x.EUTL_ID')['VERIFIED_EMISSIONS_2016']\
                                    .iloc[0], axis = 1)
    df['verified_emissions_15'] = df.apply(lambda x: ETS.query('REGISTRY_CODE==@x.EUTL_countrycode')\
                                    .query('INSTALLATION_IDENTIFIER==@x.EUTL_ID')['VERIFIED_EMISSIONS_2015']\
                                    .iloc[0], axis = 1)
    
    df['ETS_name'] = df.apply(lambda x: ETS.query('REGISTRY_CODE==@x.EUTL_countrycode')\
                                    .query('INSTALLATION_IDENTIFIER==@x.EUTL_ID')['IDENTIFIER_IN_REG']\
                                    .iloc[0], axis = 1)
    return df
generation_unit_info_matched = connect_CO2(generation_unit_info_matched, EUTL_emissions)<jupyter_output><empty_output><jupyter_text>Show first 10 entries from data<jupyter_code>generation_unit_info_matched.head(5)<jupyter_output><empty_output><jupyter_text>Carbon intensity (CI) per power plants (sites) The following function calculates the CI for each entry in the EUTL data set. Due to the fact, that the EUTL data set has a spatial resolution on site level and not on power plant level it can happen, that power plants data need to be aggregated together in the emission calculation.<jupyter_code>def calc_side_gen(df):
    # calculate the electricity generation per side.
    side_generation = df.groupby(['EUTL_countrycode','EUTL_ID']).sum()['generation_2018'].reset_index()
    
    # add side generation to PP dataset
    df['side_generation'] = df.apply(lambda x: side_generation.query('EUTL_countrycode == @x.EUTL_countrycode')\
                                                              .query('EUTL_ID == @x.EUTL_ID').iloc[0,2],axis = 1)
    
    return df
generation_unit_info_matched = calc_side_gen(generation_unit_info_matched)<jupyter_output><empty_output><jupyter_text>Distribute the reported emissions per site across the individual power plants in relation to the reported generation.<jupyter_code>generation_unit_info_matched.verified_emissions_18 = generation_unit_info_matched.verified_emissions_18.apply(float)
generation_unit_info_matched['PP_emissions'] = (generation_unit_info_matched.verified_emissions_18 / generation_unit_info_matched.side_generation) * generation_unit_info_matched.generation_2018<jupyter_output><empty_output><jupyter_text>Calculates the emission factor for each PP.<jupyter_code>generation_unit_info_matched['PP_EF'] = generation_unit_info_matched.PP_emissions / generation_unit_info_matched.generation_2018<jupyter_output><empty_output><jupyter_text>EF calculation old, EF on side level def calc_CI(df):     calculates the carbon intensity of sites. Due to non-unique labelings of the power plants across different     countries it needs care taken about it.        df['verified_emissions_18']=df['verified_emissions_18'].apply(float)        sites = pd.DataFrame(df.groupby(['EUTL_countrycode','EUTL_ID']).mean()['verified_emissions_18']\                       / df.groupby(['EUTL_countrycode','EUTL_ID']).sum()['yearly_generation_2018'])\                        .reset_index()    df['carbon_intensity'] = df.apply(lambda x: sites.query('EUTL_countrycode == @x.EUTL_countrycode')\                            .query('EUTL_ID == @x.EUTL_ID')[0].iloc[0],axis = 1)        return df<jupyter_code>#generation_unit_info_matched = calc_CI(generation_unit_info_matched)<jupyter_output><empty_output><jupyter_text>Show first 10 entries from data.Attention, the verified_emissions_* per power plant are always site-related, if there are several power plants at one site, the emissions have to be distributed among the power plants. (see. PP_emissions)<jupyter_code>generation_unit_info_matched.head()<jupyter_output><empty_output><jupyter_text>Validation of carbon intensities (CI) per power plant Check and test results and drop all results that are not plausible Which technologies are covered?<jupyter_code>generation_unit_info_matched.ProductionTypeName.unique()<jupyter_output><empty_output><jupyter_text>Missing technologies are oil and waste  covert technologies for methodetech = ['gas', 'hard_coal',       'lignite', 'other_fossil']Powerplants_emission_EU_CI = generation_unit_info_matched.query('ProductionTypeName in @tech').copy() not covert : 'Fossil Oil shale', 'Fossil Oil', 'Fossil Coal-derived gas','Fossil Peat' Since not all emission factors are correctly calculated we limited the list of calculated emissions per power plant is checked by a plausibility check.thresholds for plausibility check - until [g/kWh]: - lignite_high = 1750- lignite_low = 800-------------------- hardcoal_high = 1350- hardcoal_low = 700-------------------- gas_high = 550- gas_low = 350-------------------- other_high = 1570- other_low = 400Numbers taken from:https://www.gegenwind-saarland.de/Materialien/Energiewende/071031--VdI---CO2-Emissionen%20der%20Stromerzeugung_01.pdf<jupyter_code># emissions space per technology in [t CO2/MWh] 

emission_space = {"lignite": (1.750,.800),
                  "hard_coal": (1.350,.700),
                  "gas": (.550,.350),
                  "other_fossil": (1.570,.400)}


# Entwicklung der spezifischen Kohlendioxid-Emissionen des deutschen Strommix in den Jahren 1990 - 2018 (UBA)
# https://www.umweltbundesamt.de/sites/default/files/medien/1410/publikationen/2019-04-10_cc_10-2019_strommix_2019.pdf
#
# Unit: g CO2/kWh
#  
#lignite 1090,0
#coal 820
#gas 360
#other 1570 (Müll und Öl)
#
#
# um von t/MWh auf g/kWh zu kommen * 1000 
# function for CI check
def check_CI (row):
    if row.PP_EF > emission_space[row.ProductionTypeName][1] and row.PP_EF < emission_space[row.ProductionTypeName][0]:
        return True
    else:
        return False  
# us CI_check function for CI validation
generation_unit_info_matched['CI_validation'] = generation_unit_info_matched.apply(check_CI, axis=1)

generation_unit_info_matched.head()
len(generation_unit_info_matched)
# take only the data that are validated
generation_unit_info_validated = generation_unit_info_matched[generation_unit_info_matched.CI_validation == True]<jupyter_output><empty_output><jupyter_text>reasons why the CI is not valid:1.) Match is not correct - then one more iteration is necessary2.) Data is wrong or incomplete3.) Match is correct, but it's a highly aggregated facility such like a steel manufacturing factory, where both process-CO2 AND the emissions from the power plant running the process accounts into the number given in the dataset.4.) the facility contains more smaller power plants which are not covered by ENTSO-E (<100MW) but neverthelessproduce emissions which need to be paid for.5.) if the number of the carbon-intensity is negative, then the dataset is explained as " -1 = blank (No allocation has been made / No Emissions have been Verified)"<jupyter_code>generation_unit_info_validated.columns
len(generation_unit_info_validated)
generation_unit_info_validated.head()
fig, ax = plt.subplots(1, 1)

ax = sns.boxplot(y="ProductionTypeName", x="PP_EF", data=generation_unit_info_validated, orient="h")

ax.set_xlabel("CO2 intensity [t CO2/MWh]")
ax.set_ylabel("Production technologie")
fig.savefig(os.path.join(output_directory_path + '/_CO2_intensity_by_technology.png'))
generation_unit_info_validated.to_csv(output_directory_path + '/generation_unit_EF_bottom_up_methode.csv')<jupyter_output><empty_output><jupyter_text>Mean carbon intensity per technology and country<jupyter_code>generation_per_tech_country = generation_unit_info_validated.groupby(['EUTL_countrycode','ProductionTypeName']).sum()[['generation_2018','PP_emissions']]
EF_per_tech_country = generation_per_tech_country.PP_emissions / generation_per_tech_country.generation_2018<jupyter_output><empty_output><jupyter_text>Plot mean CI for countries and technologies<jupyter_code>plt.figure()

EF_per_tech_country.unstack(level=1).plot(kind='bar', subplots=True, ylabel='CI [t CO2/MWh]')

plt.grid(True)
EF_per_tech_country
EF_per_tech_country.to_dict()
# reset index
generation_unit_info_validated.reset_index(inplace=True)
# set new index for maping function 
generation_unit_info_EF_tech = generation_unit_info_validated.set_index(['EUTL_countrycode','ProductionTypeName'])


generation_unit_info_EF_tech['Tech_EF'] = generation_unit_info_EF_tech.index.map(EF_per_tech_country.to_dict())

#reset index
generation_unit_info_EF_tech.reset_index(inplace=True)

#set old index
generation_unit_info_EF_tech.set_index(['eic_g'],inplace=True)
generation_unit_info_EF_tech.head(5)
generation_unit_info_EF_tech.to_csv(output_directory_path + '/generation_unit_EF_bottom_up_methode.csv')<jupyter_output><empty_output><jupyter_text>Mean carbon intensity per country<jupyter_code>generation_per_country = generation_unit_info_validated.groupby(['EUTL_countrycode']).sum()[['generation_2018','PP_emissions']]
EF_per_country = generation_per_country.PP_emissions / generation_per_country.generation_2018
# reset index
generation_unit_info_EF_tech.reset_index(inplace=True)
# set new index for maping function 
generation_unit_info_EF_country = generation_unit_info_EF_tech.set_index(['MapCode'])


generation_unit_info_EF_country['Country_EF'] = generation_unit_info_EF_country.index.map(EF_per_country.to_dict())

#reset index
generation_unit_info_EF_country.reset_index(inplace=True)

#set old index
generation_unit_info_EF_country.set_index(['GenerationUnitEIC'],inplace=True)
generation_unit_info_EF_country.head(5)
generation_unit_info_EF_country.to_csv(output_directory_path + '/generation_unit_EF_bottom_up_methode.csv')<jupyter_output><empty_output><jupyter_text>Check calculated CIs (are there represent the country production) By comparing the used power plant production with the production per production type (from ENTSO-E) of each country, we can get an idea how trustful is our resulting CI.<jupyter_code>entsoe_gen_type = load_timeseries_opsd(Bootom_up_methode_input_directory_path, 'time_series_60min_multiindex.csv')<jupyter_output><empty_output><jupyter_text>Show the head of the data (the data is using a multi index for the columns)<jupyter_code>entsoe_gen_type.head()<jupyter_output><empty_output><jupyter_text>Filter the data to the countries and technologies that are covert from our method- countries- technologies<jupyter_code>entsoe_gen_type_pp = entsoe_gen_type[generation_unit_info_validated.EUTL_countrycode.unique()]
entsoe_gen_type_pp = entsoe_gen_type_pp.loc[:,(slice(None),generation_unit_info_validated.ProductionTypeName.unique())].sum()<jupyter_output><empty_output><jupyter_text>Reset the multiindex from pandas and drop not used columns<jupyter_code>entsoe_gen_type_pp = entsoe_gen_type_pp.reset_index()
entsoe_gen_type_pp = entsoe_gen_type_pp.drop(columns=['attribute', 'source', 'web', 'unit'])<jupyter_output><empty_output><jupyter_text>Rename columns<jupyter_code>entsoe_gen_type_pp = entsoe_gen_type_pp.rename(columns={0:'generation_per_type_2018_ENTSO-E','region':'EUTL_countrycode','variable':'ProductionTypeName'})<jupyter_output><empty_output><jupyter_text>Group production for each technology for the individual power plants<jupyter_code>Powerplants_gen_type_year = pd.DataFrame(generation_unit_info_validated.groupby(['EUTL_countrycode','ProductionTypeName']).sum()['generation_2018'])<jupyter_output><empty_output><jupyter_text>Reset the multiindex from pandas and drop not used columns<jupyter_code>Powerplants_gen_type_year = Powerplants_gen_type_year.reset_index()
Powerplants_gen_type_year = Powerplants_gen_type_year.rename(columns={'generation_2018':'verified_generation_per_unit_2018_ENTSO-E'})
gen_data_comparison = entsoe_gen_type_pp.merge(Powerplants_gen_type_year, on=['EUTL_countrycode','ProductionTypeName'])
gen_data_comparison.head(10)
gen_data_comparison['missing_gen'] = gen_data_comparison['generation_per_type_2018_ENTSO-E'] - gen_data_comparison['verified_generation_per_unit_2018_ENTSO-E']
gen_data_comparison['gen_covered_pro'] = round((100 / gen_data_comparison['generation_per_type_2018_ENTSO-E']) * gen_data_comparison['verified_generation_per_unit_2018_ENTSO-E'])<jupyter_output><empty_output><jupyter_text>covered generation per country and generation type<jupyter_code>ax = sns.barplot(x="EUTL_countrycode", y="gen_covered_pro", hue="ProductionTypeName", data=gen_data_comparison)
ax.plot([-1, 18], [30, 30], 'r--')<jupyter_output><empty_output><jupyter_text>Mean coverage per country<jupyter_code>fig, ax = plt.subplots(1, 1)
ax = gen_data_comparison.groupby(['EUTL_countrycode']).mean()['gen_covered_pro'].plot(kind='bar', title='Mean coverage per country', ylabel='Coverage [%]')
ax.plot([-10, 100], [30, 30], 'r--')
fig.savefig(os.path.join(output_directory_path + '/_mean_coverage_of_elec_production_per_country.png'))<jupyter_output><empty_output><jupyter_text>Mean coverage over all countries<jupyter_code>gen_data_comparison.gen_covered_pro.mean()<jupyter_output><empty_output><jupyter_text>Carbon intensity per technology and country Take only the countries and technologies where at least 30 % of the generation is connected to CO2 emissions.<jupyter_code>gen_data_comparison['Verified_for_CI_value'] = (round(gen_data_comparison['gen_covered_pro']) > 30)
gen_data_comparison.head(3)
gen_data_comparison = gen_data_comparison.set_index(['MapCode','ProductionTypeName'])
CI_PP_per_tech_country = EF_per_tech_country.reset_index().rename(columns={0:'PP_EF','countrycode':'MapCode'}).set_index(['MapCode','ProductionTypeName'])<jupyter_output><empty_output><jupyter_text>Transfer the verification of CI too the result<jupyter_code>CI_PP_per_tech_country['Verified_CI_value'] = gen_data_comparison['Verified_for_CI_value']
CI_PP_per_tech_country['EF_round'] = round(CI_PP_per_tech_country['PP_EF'],2)<jupyter_output><empty_output><jupyter_text>export the CIs to csv<jupyter_code>CI_PP_per_tech_country.to_csv(output_directory_path + '/CI_bottom_up_methode.csv')<jupyter_output><empty_output><jupyter_text>Mean carbon intensity per country<jupyter_code>EF_per_country = generation_unit_info_validated.groupby(['countrycode']).mean()['PP_EF']
EF_per_country<jupyter_output><empty_output><jupyter_text>Carbon intensity for power generation<jupyter_code>entsoe_gen_type_sum = entsoe_gen_type[Powerplants_emission_EU_CI_validated.MapCode.unique()].sum().reset_index().drop(columns=['attribute', 'source', 'web', 'unit'])
entsoe_gen_type_sum = entsoe_gen_type_sum.rename(columns={0:'generation_per_type_2018_ENTSO-E','region':'MapCode','variable':'ProductionTypeName'})
entsoe_gen_type_sum.set_index(['MapCode', 'ProductionTypeName'], inplace = True)
entsoe_gen_type_sum['emissions'] = entsoe_gen_type_sum['generation_per_type_2018_ENTSO-E'] * CI_PP_per_tech_country['CI']
entsoe_gen_type_sum = entsoe_gen_type_sum.reset_index()
CI_country = entsoe_gen_type_sum.groupby('MapCode').sum()['emissions'] / entsoe_gen_type_sum.groupby('MapCode').sum()['generation_per_type_2018_ENTSO-E']
entsoe_gen_type_sum.groupby('MapCode').sum()['generation_per_type_2018_ENTSO-E'].plot(kind='pie')
CI_country.to_csv(output_directory_path + '/CI_country_bottom_up_methode.csv')
CI_eea = pd.read_csv(output_directory_path + '/CI_top_down_methode.csv', index_col='Year')/1000
CI_eea.loc[2018]
df = pd.concat([CI_eea.loc[2018],CI_country], axis=1)
df.rename(columns={0:'Bottom up method',2018:'Top down method'}, inplace=True)

# making a nice plot:
df.plot(kind='barh', label=['Calculation','asd'])


#plt.barh(CI,label='Calculation',height=0.4, align='edge',color='b')
#plt.barh(CI_eea.loc[2018],label='EEA entry',height=-0.4,align='edge',color='y')

plt.xlabel('Carbon intensity in kg/kWh')
plt.ylabel('Countrys')
plt.title('Validation of calculated country and technology specific carbon intensities',fontweight='bold')
plt.legend()

plt.savefig(os.path.join(output_directory_path + '/_validation_of_calculation.png'))<jupyter_output><empty_output><jupyter_text>Some other checks Plausibility check for EF from power plants<jupyter_code># number of PP from matching
len(Powerplants_emission_EU_CI)
# validated PP with EF
len(Powerplants_emission_EU_CI_validated)
# covered generation by validated PP
Powerplants_emission_EU_CI_validated.yearly_generation_2018.sum()
# conv. generation in 2018
entsoe_gen_type['generation_per_type_2018_ENTSO-E'].sum()
# in % of validated PP generation
(100 / (entsoe_gen_type['generation_per_type_2018_ENTSO-E'].sum())) * Powerplants_emission_EU_CI_validated.yearly_generation_2018.sum()
pl_match = unit_matching_EU[unit_matching_EU.countrycode == 'PL']
pl_match[pl_match.EUTL_ID == 1]
Powerplants_emission_EU_CI[Powerplants_emission_EU_CI.PowerSystemResourceName.str.contains("Bełchatów")]
generation_per_unit
generation_unit_info[generation_unit_info.PowerSystemResourceName.str.contains("Bełchatów")]
PL = EUTL_emissions[EUTL_emissions.REGISTRY_CODE == 'PL']
PL = PL[PL.VERIFIED_EMISSIONS_2018 != 'Excluded']
PL.VERIFIED_EMISSIONS_2018 = PL.VERIFIED_EMISSIONS_2018.astype(float)
PL.sort_values(by='VERIFIED_EMISSIONS_2018', ascending=False)
coverage_val = len(Powerplants_emission_EU_CI) / len(generation_unit_info.query('ProductionTypeName in @tech'))
print("After validation we were able to allocate " + str(coverage_val*100) + " % of the power plants an emission value.")
Powerplants_emission_EU_CI.yearly_generation_2018.sum()

coverage = len(Powerplants_emission_EU_CI) / len(generation_unit_info.query('ProductionTypeName in @tech'))
print("We were able to allocate " + str(coverage*100) + " % of the power plants an emission value.")
# Plot CI per country for the technologies (drei Plots)
Powerplants_emission_EU_CI.columns
# Teste wie viel Energie decken wir mit dieser Methode ab?
# Dafür welche Daten verwenden (Problem Gross vs Net)
# Gibt es von Eurostat eine Net Version (ich glaube ja pro Monat)
# Dann Plot pro Land wie viel Dekcen wir ab (0 bis 1)


# wie könnten wir ein representatives sample bilder. Das soll uns dann die CI ausrechnen

# Wie machen wir den Vergleich zur anderen Methode?
# Jetzt können wir dann CI pro technologie bestimmen.<jupyter_output><empty_output><jupyter_text>vergleich unsere Jährliche Erzeugung mit denen, die wir von Ember bekommen haben (in 12, in Gfuellner Code) . Die unserscheidet sich. Was macht Ember noch mit den Daten?<jupyter_code>#generation_unit_info[generation_unit_info.PowerSystemResourceName == 'Scheldelaan Exxonmobil']
pd.merge(generation_unit_info, unit_matching_EU, on='PowerSystemResourceName', how='left')<jupyter_output><empty_output><jupyter_text>Some check up numbers<jupyter_code>tech = ['Fossil Gas', 'Fossil Hard coal',
       'Fossil Brown coal/Lignite', 'Other']

print('number of pp in data set that has emissions that we want to match')
print(len(generation_unit_info.query('ProductionTypeName in @tech')))

print('generation in MWh in 2018 from the unis')
print(generation_unit_info.query('ProductionTypeName in @tech').yearly_generation_2018.sum())


(100/(generation_unit_info.query('ProductionTypeName in @tech').yearly_generation_2018.sum()))* Powerplants_emission_EU_CI.yearly_generation_2018.sum()
generation_unit_info_matched[generation_unit_info_matched.ProductionTypeName == 'Other'].sort_values(by='yearly_generation_2018', ascending=False)
generation_unit_info.ProductionTypeName.unique()
# überprüfen, ob wir noch weitere technologien mit aufnehmen sollen 
# ich denke IT 'other' sollte man sich anschauen da haben wir auch schon einige drin.
tech = ['Fossil Oil shale', 'Fossil Oil','Fossil Coal-derived gas',
       'Fossil Peat', 'Other', 'Waste']
generation_unit_info.query('ProductionTypeName in @tech').sort_values(by='yearly_generation_2018', ascending=False).head(50)<jupyter_output><empty_output><jupyter_text>--- Following function and data are only for check up the results and are not needed (test data)  ---  1) Entso unit generation data from EMBERThe Entso e data has lot of gaps and need to be preprocessing (is done partly and the result is different than the one from Ember)File - > entsoe_unit_gen_data.csv<jupyter_code>def load_unit_data_ember(path, fn):
    """
    Entso e unit generation from ember
        
    Parameters
    ----------
    path: str
        path to data
    fn : str
        filename
        
    """
    
    logging.info(f'Loading data from {fn}')
    
    df = pd.read_csv(os.path.join(path, fn), sep = ';', header = 0, encoding = 'unicode_escape')

    return df

# load the ember data
generation_per_unit_ember = load_unit_data_ember(Bootom_up_methode_input_directory_path, 'entsoe_unit_gen_data.csv') 

# filter to the year 2018
generation_per_unit_ember = generation_per_unit_ember[generation_per_unit_ember.Year == 2018]

# group and sum monthly data too yearly data
generation_per_unit_ember = generation_per_unit_ember.groupby(by = 'PowerSystemResourceName').sum()['Monthly Generation (GWh)']
# for the ember datat the country code in germany needs to be converted 
def conv(x):
    # converts german grid operator areas string in german countrycode of ETS
    areas =['DE_TenneT_GER','DE_TransnetBW','DE_Amprion','DE_50HzT']
    if x in areas:
        x = 'DE' 
    return x
#not needed anymore
#generation_per_unit_ember = generation_per_unit_ember.groupby(by = 'PowerSystemResourceName').agg({'Monthly Generation (GWh)': 'sum',
                                                                       'MapCode':'first',
                                                                       'ProductionTypeName':'first',
                                                                       
                                                                         'Year': 'first',
                                                                         'Installed Capacity (MW)': 'first'
                                                                      })
#merged_ember = pd.merge(generation_per_unit_ember, unit_matching_EU_2, on='PowerSystemResourceName', how='inner')<jupyter_output><empty_output><jupyter_text>CO2 SignalsIn this section, we calculate CO2 signals based on CO2 emissions from power plant.   On power plant level Using the above calculated EF per power plant to design an hourly CO2 signal per power plant.<jupyter_code>hourly_CO2_emission_per_unit = generation_per_unit.multiply((Powerplants_emission_EU_CI_validated.carbon_intensity), axis='columns')<jupyter_output><empty_output><jupyter_text>Drop all power plants without an EF value.<jupyter_code>hourly_CO2_emission_per_unit.dropna(axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>Filter Data to Germany availabe power plants in Germany<jupyter_code>DE_PP_ID = Powerplants_emission_EU_CI_validated[Powerplants_emission_EU_CI_validated.MapCode == 'DE'].index.to_list()<jupyter_output><empty_output><jupyter_text>German Dataset<jupyter_code>hourly_CO2_emission_per_unit_DE = hourly_CO2_emission_per_unit[DE_PP_ID]
hourly_CO2_emission_per_unit_DE.sum(axis=1).plot()
(hourly_CO2_emission_per_unit_DE.sum(axis=1) / generation_per_unit[DE_PP_ID].sum(axis=1)).describe()<jupyter_output><empty_output><jupyter_text>some tests<jupyter_code>#CO2_per_unit[Powerplants_emission_EU_CI_validated.index].sum().loc[['48W00000EGGPS-1X','48W00000EGGPS-4R','48W00000EGGPS-3T','48W00000EGGPS-2V']].sum()
#Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.countrycode == 'GB') & (Powerplants_emission_EU_CI_validated.EUTL_ID == 169)]<jupyter_output><empty_output><jupyter_text>On tech level in Germany<jupyter_code>Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')].head()
def calc_CI_tech_country_new(df):
    # calculates the average carbon intensity per technology in a country for plausible results defined above.
    # df is a dataframe which contains all sites with plausible carbon intensities of whole Europe.
    CI_tech_country = df.groupby(['countrycode','EUTL_ID']).agg({'verified_emissions_18':'mean',
                                                                 'ProductionTypeName':'first'}).groupby(['ProductionTypeName']).sum()['verified_emissions_18']\
                                     /\
                      df.groupby(['countrycode','EUTL_ID']).agg({'yearly_generation_2018':'sum',
                                                                 'ProductionTypeName':'first'}).groupby(['ProductionTypeName']).sum()['yearly_generation_2018']

    return CI_tech_country
EF_tech_dic = calc_CI_tech_country_new(Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')]).to_dict()
generation_unit_info['EF'] = generation_unit_info.ProductionTypeName.map(EF_tech_dic)
generation_unit_info.head()
dic = generation_unit_info.EF.to_dict()
dic
generation_per_unit_DE = generation_per_unit[Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')].index.to_list()]
hourly_CO2_emission_per_unit_2 = generation_per_unit_DE.assign(**dic).mul(generation_per_unit_DE)
hourly_CO2_emission_per_unit_DE2 = hourly_CO2_emission_per_unit_2[DE_PP_ID]
(hourly_CO2_emission_per_unit_DE2.sum(axis=1) - hourly_CO2_emission_per_unit_DE.sum(axis=1)).plot()<jupyter_output><empty_output><jupyter_text>On country level in Germany<jupyter_code>def calc_CI_country_new(df):
    # calculates the average carbon intensity per technology in a country for plausible results defined above.
    # df is a dataframe which contains all sites with plausible carbon intensities of whole Europe.
    CI_country = df.groupby(['countrycode','EUTL_ID']).agg({'verified_emissions_18':'mean',
                                                                 'ProductionTypeName':'first'}).groupby(['countrycode']).sum()['verified_emissions_18']\
                                     /\
                      df.groupby(['countrycode','EUTL_ID']).agg({'yearly_generation_2018':'sum',
                                                                 'ProductionTypeName':'first'}).groupby(['countrycode']).sum()['yearly_generation_2018']

    return CI_country
EF_dic = calc_CI_country_new(Powerplants_emission_EU_CI_validated).to_dict()
EF_dic
generation_unit_info['EF'] = generation_unit_info.MapCode.map(EF_dic)
dic = generation_unit_info.EF.to_dict()
generation_per_unit_DE = generation_per_unit[Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')].index.to_list()]
hourly_CO2_emission_per_unit_3 = generation_per_unit_DE.assign(**dic).mul(generation_per_unit_DE)
hourly_CO2_emission_per_unit_DE3 = hourly_CO2_emission_per_unit_3[DE_PP_ID]
(hourly_CO2_emission_per_unit_DE3.sum(axis=1) - hourly_CO2_emission_per_unit_DE.sum(axis=1)).plot()
EF = pd.DataFrame()
EF['PP'] = hourly_CO2_emission_per_unit_DE.sum(axis=1)
EF['tech'] = hourly_CO2_emission_per_unit_DE2.sum(axis=1)
EF['country'] = hourly_CO2_emission_per_unit_DE3.sum(axis=1)
EF.plot.scatter(x='tech', y='country')<jupyter_output><empty_output><jupyter_text>hier ist es jetzt natürlich spannend auch EE Anlagen mit in die Rechnung mit rein zu nehmen, da sie einen EF von 0 haben, was sich bei tech und PP gut zeigen würde, bei country aber nicht.<jupyter_code>EF.describe()
sns.boxplot(y="ProductionTypeName", x="carbon_intensity", data=EF, orient="h")<jupyter_output><empty_output><jupyter_text>CO2 SignalsIn this section, we calculate CO2 signals based on CO2 emissions from power plant.   On power plant level Using the above calculated EF per power plant to design an hourly CO2 signal per power plant.<jupyter_code>hourly_CO2_emission_per_unit = generation_per_unit.multiply((Powerplants_emission_EU_CI_validated.carbon_intensity), axis='columns')<jupyter_output><empty_output><jupyter_text>Drop all power plants without an EF value.<jupyter_code>hourly_CO2_emission_per_unit.dropna(axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>Filter Data to Germany availabe power plants in Germany<jupyter_code>DE_PP_ID = Powerplants_emission_EU_CI_validated[Powerplants_emission_EU_CI_validated.MapCode == 'DE'].index.to_list()<jupyter_output><empty_output><jupyter_text>German Dataset<jupyter_code>hourly_CO2_emission_per_unit_DE = hourly_CO2_emission_per_unit[DE_PP_ID]
hourly_CO2_emission_per_unit_DE.sum(axis=1).plot()
(hourly_CO2_emission_per_unit_DE.sum(axis=1) / generation_per_unit[DE_PP_ID].sum(axis=1)).describe()<jupyter_output><empty_output><jupyter_text>some tests<jupyter_code>#CO2_per_unit[Powerplants_emission_EU_CI_validated.index].sum().loc[['48W00000EGGPS-1X','48W00000EGGPS-4R','48W00000EGGPS-3T','48W00000EGGPS-2V']].sum()
#Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.countrycode == 'GB') & (Powerplants_emission_EU_CI_validated.EUTL_ID == 169)]<jupyter_output><empty_output><jupyter_text>On tech level in Germany<jupyter_code>Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')].head()
def calc_CI_tech_country_new(df):
    # calculates the average carbon intensity per technology in a country for plausible results defined above.
    # df is a dataframe which contains all sites with plausible carbon intensities of whole Europe.
    CI_tech_country = df.groupby(['countrycode','EUTL_ID']).agg({'verified_emissions_18':'mean',
                                                                 'ProductionTypeName':'first'}).groupby(['ProductionTypeName']).sum()['verified_emissions_18']\
                                     /\
                      df.groupby(['countrycode','EUTL_ID']).agg({'yearly_generation_2018':'sum',
                                                                 'ProductionTypeName':'first'}).groupby(['ProductionTypeName']).sum()['yearly_generation_2018']

    return CI_tech_country
EF_tech_dic = calc_CI_tech_country_new(Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')]).to_dict()
generation_unit_info['EF'] = generation_unit_info.ProductionTypeName.map(EF_tech_dic)
generation_unit_info.head()
dic = generation_unit_info.EF.to_dict()
dic
generation_per_unit_DE = generation_per_unit[Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')].index.to_list()]
hourly_CO2_emission_per_unit_2 = generation_per_unit_DE.assign(**dic).mul(generation_per_unit_DE)
hourly_CO2_emission_per_unit_DE2 = hourly_CO2_emission_per_unit_2[DE_PP_ID]
(hourly_CO2_emission_per_unit_DE2.sum(axis=1) - hourly_CO2_emission_per_unit_DE.sum(axis=1)).plot()<jupyter_output><empty_output><jupyter_text>On country level in Germany<jupyter_code>def calc_CI_country_new(df):
    # calculates the average carbon intensity per technology in a country for plausible results defined above.
    # df is a dataframe which contains all sites with plausible carbon intensities of whole Europe.
    CI_country = df.groupby(['countrycode','EUTL_ID']).agg({'verified_emissions_18':'mean',
                                                                 'ProductionTypeName':'first'}).groupby(['countrycode']).sum()['verified_emissions_18']\
                                     /\
                      df.groupby(['countrycode','EUTL_ID']).agg({'yearly_generation_2018':'sum',
                                                                 'ProductionTypeName':'first'}).groupby(['countrycode']).sum()['yearly_generation_2018']

    return CI_country
EF_dic = calc_CI_country_new(Powerplants_emission_EU_CI_validated).to_dict()
EF_dic
generation_unit_info['EF'] = generation_unit_info.MapCode.map(EF_dic)
dic = generation_unit_info.EF.to_dict()
generation_per_unit_DE = generation_per_unit[Powerplants_emission_EU_CI_validated[(Powerplants_emission_EU_CI_validated.MapCode == 'DE')].index.to_list()]
hourly_CO2_emission_per_unit_3 = generation_per_unit_DE.assign(**dic).mul(generation_per_unit_DE)
hourly_CO2_emission_per_unit_DE3 = hourly_CO2_emission_per_unit_3[DE_PP_ID]
(hourly_CO2_emission_per_unit_DE3.sum(axis=1) - hourly_CO2_emission_per_unit_DE.sum(axis=1)).plot()
EF = pd.DataFrame()
EF['PP'] = hourly_CO2_emission_per_unit_DE.sum(axis=1)
EF['tech'] = hourly_CO2_emission_per_unit_DE2.sum(axis=1)
EF['country'] = hourly_CO2_emission_per_unit_DE3.sum(axis=1)
EF.plot.scatter(x='tech', y='country')<jupyter_output><empty_output><jupyter_text>hier ist es jetzt natürlich spannend auch EE Anlagen mit in die Rechnung mit rein zu nehmen, da sie einen EF von 0 haben, was sich bei tech und PP gut zeigen würde, bei country aber nicht.<jupyter_code>EF.describe()
sns.boxplot(y="ProductionTypeName", x="carbon_intensity", data=EF, orient="h")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10439037966463344, 'fraction_numerical': 0.016768328037664957, 'mean_word_length': 5.221035925551869, 'pattern_counts': {'":': 4, '<': 234, '<?xml version=': 0, '>': 247, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 49, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '58', 'dir': 'jupyter-structured-clean-dedup', 'id': '180933', 'license': 'MIT', 'n_tokens_mistral': 14936, 'n_tokens_neox': 13759, 'n_words': 2655, 'path': '.ipynb_checkpoints/CO2 calculation bottom up method-checkpoint.ipynb', 'repo_name': 'INATECH-CIG/CO2_emissions_factors', 'size': '43117'}