data.datasets.emobility.motorized_individual_travel.ev_allocation.calc_evs_per_grid_district() - Code Metrics - Inspection of "Merge pull request #1138 from openego/documentatio..." - openego/eGon-data - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — dev ( 85d654...312d71 )

unknown

created 2023-08-31 16:57 UTC

calc_evs_per_grid_district() B

↳ Parent: data.datasets.emobility.motorized_individual_travel.ev_allocation

Complexity

Conditions

Size

Total Lines	130
Code Lines	77

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	77
dl	0
loc	130
rs	7.7527
c	0
b	0
f	0
cc	4
nop	1

How to fix Long Method

"""
* Calculate number of electric vehicles and allocate on different spatial
  levels: :py:func:`allocate_evs_numbers`
* Allocate specific EVs to MV grid districts:
  :py:func:`allocate_evs_to_grid_districts`

"""

from itertools import permutations

from sqlalchemy.sql import func
import numpy as np
import pandas as pd

from egon.data import db
from egon.data.datasets.emobility.motorized_individual_travel.db_classes import (
    EgonEvCountMunicipality,
    EgonEvCountMvGridDistrict,
    EgonEvCountRegistrationDistrict,
    EgonEvMvGridDistrict,
    EgonEvPool,
)
from egon.data.datasets.emobility.motorized_individual_travel.helpers import (
    COLUMNS_KBA,
    CONFIG_EV,
    TESTMODE_OFF,
    read_kba_data,
    read_rs7_data,
)
from egon.data.datasets.emobility.motorized_individual_travel.tests import (
    validate_electric_vehicles_numbers,
)
from egon.data.datasets.mv_grid_districts import MvGridDistricts
from egon.data.datasets.scenario_parameters import get_sector_parameters
from egon.data.datasets.zensus_mv_grid_districts import MapZensusGridDistricts
from egon.data.datasets.zensus_vg250 import (
    DestatisZensusPopulationPerHaInsideGermany,
    MapZensusVg250,
    Vg250Gem,
    Vg250GemPopulation,
)
import egon.data.config

RANDOM_SEED = egon.data.config.settings()["egon-data"]["--random-seed"]


def fix_missing_ags_municipality_regiostar(muns, rs7_data):
    """Check if all AGS of municipality dataset are included in RegioStaR7
    dataset and vice versa.

    As of Dec 2021, some municipalities are not included int the RegioStaR7
    dataset. This is mostly caused by incorporations of a municipality by
    another municipality. This is fixed by assigning a RS7 id from another
    municipality with similar AGS (most likely a neighboured one).

    Missing entries in the municipality dataset is printed but not fixed
    as it doesn't result in bad data. Nevertheless, consider to update the
    municipality/VG250 dataset.

    Parameters
    ----------
    muns : pandas.DataFrame
        Municipality data
    rs7_data : pandas.DataFrame
        RegioStaR7 data

    Returns
    -------
    pandas.DataFrame
        Fixed RegioStaR7 data
    """

    if len(muns.ags.unique()) != len(rs7_data.ags):
        print(
            "==========> Number of AGS differ between VG250 and RS7, "
            "trying to fix this..."
        )

        # Get AGS differences
        ags_datasets = {"RS7": rs7_data.ags, "VG250": muns.ags}
        ags_datasets_missing = {k: [] for k in ags_datasets.keys()}
        perm = permutations(ags_datasets.items())
        for (name1, ags1), (name2, ags2) in perm:
            print(f"  Checking if all AGS of {name1} are in {name2}...")
            missing = [_ for _ in ags1 if _ not in ags2.to_list()]
            if len(missing) > 0:
                ags_datasets_missing[name2] = missing
                print(f"    AGS in {name1} but not in {name2}: ", missing)
            else:
                print("    OK")

        print("")

        # Try to fix
        for name, missing in ags_datasets_missing.items():
            if len(missing) > 0:
                # RS7 entries missing: use RS7 number from mun with similar AGS
                if name == "RS7":
                    for ags in missing:
                        similar_entry = rs7_data[
                            round((rs7_data.ags.div(10))) == round(ags / 10)
                        ].iloc[0]
                        if len(similar_entry) > 0:
                            print(
                                f"Adding dataset from VG250 to RS7 "
                                f"based upon AGS {ags}."
                            )
                            similar_entry.ags = ags
                            rs7_data = rs7_data.append(similar_entry)
                        print("Consider to update RS7.")
                # VG250 entries missing:
                elif name == "VG250":
                    print(
                        "Cannot guess VG250 entries. This error does not "
                        "result in bad data but consider to update VG250."
                    )

        if len(muns.ags.unique()) != len(rs7_data.ags):
            print("==========> AGS could not be fixed!")
        else:
            print("==========> AGS were fixed!")

    return rs7_data


def calc_evs_per_reg_district(scenario_variation_parameters, kba_data):
    """Calculate EVs per registration district

    Parameters
    ----------
    scenario_variation_parameters : dict
        Parameters of scenario variation
    kba_data : pandas.DataFrame
        Vehicle registration data for registration district

    Returns
    -------
    pandas.DataFrame
        EVs per registration district
    """

    scenario_variation_parameters["mini_share"] = (
        scenario_variation_parameters["bev_mini_share"]
        + scenario_variation_parameters["phev_mini_share"]
    )
    scenario_variation_parameters["medium_share"] = (
        scenario_variation_parameters["bev_medium_share"]
        + scenario_variation_parameters["phev_medium_share"]
    )
    scenario_variation_parameters["luxury_share"] = (
        scenario_variation_parameters["bev_luxury_share"]
        + scenario_variation_parameters["phev_luxury_share"]
    )

    factor_dict = dict()
    factor_dict["mini_factor"] = (
        scenario_variation_parameters["mini_share"]
        * scenario_variation_parameters["ev_count"]
        / kba_data.mini.sum()
    )
    factor_dict["medium_factor"] = (
        scenario_variation_parameters["medium_share"]
        * scenario_variation_parameters["ev_count"]
        / kba_data.medium.sum()
    )
    factor_dict["luxury_factor"] = (
        scenario_variation_parameters["luxury_share"]
        * scenario_variation_parameters["ev_count"]
        / kba_data.luxury.sum()
    )

    # Define shares and factors
    ev_data = kba_data.copy()

    for tech, params in CONFIG_EV.items():
        ev_data[tech] = (
            (
                kba_data[params["column"]]
                * factor_dict[params["factor"]]
                * scenario_variation_parameters[params["tech_share"]]
                / scenario_variation_parameters[params["share"]]
            )
            .round()
            .astype("int")
        )

    ev_data.drop(
        columns=[_ for _ in COLUMNS_KBA if _ != "reg_district"], inplace=True
    )

    return ev_data


def calc_evs_per_municipality(ev_data, rs7_data):
    """Calculate EVs per municipality

    Parameters
    ----------
    ev_data : pandas.DataFrame
        EVs per regstration district
    rs7_data : pandas.DataFrame
        RegioStaR7 data
    """
    with db.session_scope() as session:
        query = session.query(
            Vg250GemPopulation.ags_0.label("ags"),
            Vg250GemPopulation.gen,
            Vg250GemPopulation.population_total.label("pop"),
        )

    muns = pd.read_sql(
        query.statement, query.session.bind, index_col=None
    ).astype({"ags": "int64"})

    muns["ags_district"] = (
        muns.ags.multiply(1 / 1000).apply(np.floor).astype("int")
    )

    # Manual fix of Trier-Saarburg: Introduce new `ags_reg_district`
    # for correct allocation of mun to registration district
    # (Zulassungsbezirk), see above for background.
    muns["ags_reg_district"] = muns["ags_district"]
    muns.loc[muns["ags_reg_district"] == 7235, "ags_reg_district"] = 7211

    # Remove multiple municipality entries (due to 'gf' in VG250)
    # by summing up population
    muns = (
        muns[["ags", "gen", "ags_reg_district", "pop"]]
        .groupby(["ags", "gen", "ags_reg_district"])
        .sum()
        .reset_index()
    )

    # Add population of registration district
    pop_per_reg_district = (
        muns[["ags_reg_district", "pop"]]
        .groupby("ags_reg_district")
        .sum()
        .rename(columns={"pop": "pop_district"})
        .reset_index()
    )

    # Fix missing ags in mun data if not in testmode
    if TESTMODE_OFF:
        rs7_data = fix_missing_ags_municipality_regiostar(muns, rs7_data)

    # Merge municipality, EV data and pop per district
    ev_data_muns = muns.merge(ev_data, on="ags_reg_district").merge(
        pop_per_reg_district, on="ags_reg_district"
    )

    # Disaggregate EV numbers to municipality
    for tech in ev_data[CONFIG_EV.keys()]:
        ev_data_muns[tech] = round(
            ev_data_muns[tech]
            * ev_data_muns["pop"]
            / ev_data_muns["pop_district"]
        ).astype("int")

    # Filter columns
    cols = ["ags"]
    cols.extend(CONFIG_EV.keys())
    ev_data_muns = ev_data_muns[cols]

    # Merge RS7 data
    ev_data_muns = ev_data_muns.merge(rs7_data[["ags", "rs7_id"]], on="ags")

    return ev_data_muns


def calc_evs_per_grid_district(ev_data_muns):
    """Calculate EVs per grid district by using population weighting

    Parameters
    ----------
    ev_data_muns : pandas.DataFrame
        EV data for municipalities

    Returns
    -------
    pandas.DataFrame
        EV data for grid districts
    """

    # Read MVGDs with intersecting muns and aggregate pop for each
    # municipality part
    with db.session_scope() as session:
        query_pop_per_mvgd = (
            session.query(
                MvGridDistricts.bus_id,
                Vg250Gem.ags,
                func.sum(
                    DestatisZensusPopulationPerHaInsideGermany.population
                ).label("pop"),
            )
            .select_from(MapZensusGridDistricts)
            .join(
                MvGridDistricts,
                MapZensusGridDistricts.bus_id == MvGridDistricts.bus_id,
            )
            .join(
                DestatisZensusPopulationPerHaInsideGermany,
                MapZensusGridDistricts.zensus_population_id
                == DestatisZensusPopulationPerHaInsideGermany.id,
            )
            .join(
                MapZensusVg250,
                MapZensusGridDistricts.zensus_population_id
                == MapZensusVg250.zensus_population_id,
            )
            .join(
                Vg250Gem, MapZensusVg250.vg250_municipality_id == Vg250Gem.id
            )
            .group_by(MvGridDistricts.bus_id, Vg250Gem.ags)
            .order_by(Vg250Gem.ags)
        )

    mvgd_pop_per_mun = pd.read_sql(
        query_pop_per_mvgd.statement,
        query_pop_per_mvgd.session.bind,
        index_col=None,
    ).astype({"bus_id": "int64", "pop": "int64", "ags": "int64"})

    # Calc population share of each municipality in MVGD
    mvgd_pop_per_mun_in_mvgd = mvgd_pop_per_mun.groupby(["bus_id", "ags"]).agg(
        {"pop": "sum"}
    )

    # Calc relative and absolute population shares:
    # * pop_mun_in_mvgd: pop share of mun which intersects with MVGD
    # * pop_share_mun_in_mvgd: relative pop share of mun which
    #   intersects with MVGD
    # * pop_mun_total: total pop of mun
    # * pop_mun_in_mvgd_of_mun_total: relative pop share of mun which
    #   intersects with MVGD in relation to total pop of mun
    mvgd_pop_per_mun_in_mvgd = (
        mvgd_pop_per_mun_in_mvgd.groupby(level=0)
        .apply(lambda x: x / float(x.sum()))
        .reset_index()
        .rename(columns={"pop": "pop_share_mun_in_mvgd"})
        .merge(
            mvgd_pop_per_mun_in_mvgd.reset_index(),
            on=["bus_id", "ags"],
            how="left",
        )
        .rename(columns={"pop": "pop_mun_in_mvgd"})
        .merge(
            mvgd_pop_per_mun[["ags", "pop"]]
            .groupby("ags")
            .agg({"pop": "sum"}),
            on="ags",
            how="left",
        )
        .rename(columns={"pop": "pop_mun_total"})
    )
    mvgd_pop_per_mun_in_mvgd["pop_mun_in_mvgd_of_mun_total"] = (
        mvgd_pop_per_mun_in_mvgd["pop_mun_in_mvgd"]
        / mvgd_pop_per_mun_in_mvgd["pop_mun_total"]
    )

    # Merge EV data
    ev_data_mvgds = mvgd_pop_per_mun_in_mvgd.merge(
        ev_data_muns, on="ags", how="left"
    ).sort_values(["bus_id", "ags"])

    # Calc EVs per MVGD by using EV from mun and share of mun's pop
    # that is located within MVGD
    for tech in ev_data_mvgds[CONFIG_EV.keys()]:
        ev_data_mvgds[tech] = (
            round(
                ev_data_mvgds[tech]
                * ev_data_mvgds["pop_mun_in_mvgd_of_mun_total"]
            )
            .fillna(0)
            .astype("int")
        )

    # Set RS7 id for MVGD by using the RS7 id from the mun with the
    # highest share in population
    rs7_data_mvgds = (
        ev_data_mvgds[["bus_id", "pop_mun_in_mvgd", "rs7_id"]]
        .groupby(["bus_id", "rs7_id"])
        .sum()
        .sort_values(
            ["bus_id", "pop_mun_in_mvgd"], ascending=False, na_position="last"
        )
        .reset_index()
        .drop_duplicates("bus_id", keep="first")[["bus_id", "rs7_id"]]
    )

    # Join RS7 id and select columns
    columns = ["bus_id"] + [_ for _ in CONFIG_EV.keys()]
    ev_data_mvgds = (
        ev_data_mvgds[columns]
        .groupby("bus_id")
        .agg("sum")
        .merge(rs7_data_mvgds, on="bus_id", how="left")
    )

    return ev_data_mvgds


def allocate_evs_numbers():
    """Allocate electric vehicles to different spatial levels.

    Accocation uses today's vehicles registration data per registration
    district from KBA and scales scenario's EV targets (BEV and PHEV)
    linearly using population. Furthermore, a RegioStaR7 code (BMVI) is
    assigned.

    Levels:
    * districts of registration
    * municipalities
    * grid districts

    Parameters
    ----------

    Returns
    -------

    """
    # Import
    kba_data = read_kba_data()
    rs7_data = read_rs7_data()

    for scenario_name in ["eGon2035", "eGon100RE"]:
        # Load scenario params
        scenario_parameters = get_sector_parameters(
            "mobility", scenario=scenario_name
        )["motorized_individual_travel"]

        print(f"SCENARIO: {scenario_name}")

        # Go through scenario variations
        for (
            scenario_variation_name,
            scenario_variation_parameters,
        ) in scenario_parameters.items():

            print(f"  SCENARIO VARIATION: {scenario_variation_name}")

            # Get EV target
            ev_target = scenario_variation_parameters["ev_count"]

            #####################################
            # EV data per registration district #
            #####################################
            print("Calculate EV numbers for registration districts...")
            ev_data = calc_evs_per_reg_district(
                scenario_variation_parameters, kba_data
            )
            # Check EV results if not in testmode
            if TESTMODE_OFF:
                validate_electric_vehicles_numbers(
                    "EV count in registration districts", ev_data, ev_target
                )
            # Add scenario columns and write to DB
            ev_data["scenario"] = scenario_name
            ev_data["scenario_variation"] = scenario_variation_name
            ev_data.sort_values(
                ["scenario", "scenario_variation", "ags_reg_district"],
                inplace=True,
            )
            ev_data.to_sql(
                name=EgonEvCountRegistrationDistrict.__table__.name,
                schema=EgonEvCountRegistrationDistrict.__table__.schema,
                con=db.engine(),
                if_exists="append",
                index=False,
            )

            #####################################
            #     EV data per municipality      #
            #####################################
            print("Calculate EV numbers for municipalities...")
            ev_data_muns = calc_evs_per_municipality(ev_data, rs7_data)
            # Check EV results if not in testmode
            if TESTMODE_OFF:
                validate_electric_vehicles_numbers(
                    "EV count in municipalities", ev_data_muns, ev_target
                )
            # Add scenario columns and write to DB
            ev_data_muns["scenario"] = scenario_name
            ev_data_muns["scenario_variation"] = scenario_variation_name
            ev_data_muns.sort_values(
                ["scenario", "scenario_variation", "ags"], inplace=True
            )
            ev_data_muns.to_sql(
                name=EgonEvCountMunicipality.__table__.name,
                schema=EgonEvCountMunicipality.__table__.schema,
                con=db.engine(),
                if_exists="append",
                index=False,
            )

            #####################################
            #     EV data per grid district     #
            #####################################
            print("Calculate EV numbers for grid districts...")
            ev_data_mvgds = calc_evs_per_grid_district(ev_data_muns)
            # Check EV results if not in testmode
            if TESTMODE_OFF:
                validate_electric_vehicles_numbers(
                    "EV count in grid districts", ev_data_mvgds, ev_target
                )
            # Add scenario columns and write to DB
            ev_data_mvgds["scenario"] = scenario_name
            ev_data_mvgds["scenario_variation"] = scenario_variation_name
            ev_data_mvgds.sort_values(
                ["scenario", "scenario_variation", "bus_id"], inplace=True
            )
            ev_data_mvgds.to_sql(
                name=EgonEvCountMvGridDistrict.__table__.name,
                schema=EgonEvCountMvGridDistrict.__table__.schema,
                con=db.engine(),
                if_exists="append",
                index=False,
            )


def allocate_evs_to_grid_districts():
    """Allocate EVs to MV grid districts for all scenarios and scenario
    variations.

    Each grid district in
    :class:`egon.data.datasets.mv_grid_districts.MvGridDistricts`
    is assigned a list of electric vehicles from the EV pool in
    :class:`EgonEvPool` based on the RegioStar7 region and the
    counts per EV type in :class:`EgonEvCountMvGridDistrict`.
    Results are written to :class:`EgonEvMvGridDistrict`.
    """

    def get_random_evs(row):
        """Get random EV sample for EV type and RS7 region"""
        return (
            ev_pool.loc[

                (ev_pool.rs7_id == row.rs7_id)
                & (ev_pool["type"] == row["type"])
            ]
            .sample(row["count"], replace=True)
            .ev_id.to_list()
        )

    for scenario_name in ["eGon2035", "eGon100RE"]:
        print(f"SCENARIO: {scenario_name}")

        # Load EVs per grid district
        print("Loading EV counts for grid districts...")
        with db.session_scope() as session:
            query = session.query(EgonEvCountMvGridDistrict).filter(
                EgonEvCountMvGridDistrict.scenario == scenario_name
            )
        ev_per_mvgd = pd.read_sql(
            query.statement, query.session.bind, index_col=None
        )

        # Convert EV types' wide to long format
        ev_per_mvgd = pd.melt(
            ev_per_mvgd,
            id_vars=["scenario", "scenario_variation", "bus_id", "rs7_id"],
            value_vars=CONFIG_EV.keys(),
            var_name="type",
            value_name="count",
        )

        # Load EV pool
        print("  Loading EV pool...")
        with db.session_scope() as session:
            query = session.query(EgonEvPool).filter(
                EgonEvPool.scenario == scenario_name
            )
        ev_pool = pd.read_sql(
            query.statement,
            query.session.bind,
            index_col=None,
        )

        # Draw EVs randomly for each grid district from pool
        print("  Draw EVs from pool for grid districts...")
        np.random.seed(RANDOM_SEED)
        ev_per_mvgd["egon_ev_pool_ev_id"] = ev_per_mvgd.apply(
            get_random_evs, axis=1
        )
        ev_per_mvgd.drop(columns=["rs7_id", "type", "count"], inplace=True)

        # EV lists to rows
        ev_per_mvgd = ev_per_mvgd.explode("egon_ev_pool_ev_id")

        # Check for empty entries
        empty_ev_entries = ev_per_mvgd.egon_ev_pool_ev_id.isna().sum()
        if empty_ev_entries > 0:
            print("====================================================")
            print(
                f"WARNING: Found {empty_ev_entries} empty entries "
                f"and will remove it:"
            )
            print(ev_per_mvgd[ev_per_mvgd.egon_ev_pool_ev_id.isna()])
            ev_per_mvgd = ev_per_mvgd[~ev_per_mvgd.egon_ev_pool_ev_id.isna()]
            print("====================================================")

        # Write trips to DB
        print("  Writing allocated data to DB...")
        ev_per_mvgd.to_sql(
            name=EgonEvMvGridDistrict.__table__.name,
            schema=EgonEvMvGridDistrict.__table__.schema,
            con=db.engine(),
            if_exists="append",
            index=False,
            method="multi",
            chunksize=10000,
        )

        # Check EV result sums for all scenario variations if not in testmode
        if TESTMODE_OFF:
            print("  Validating results...")
            ev_per_mvgd_counts_per_scn = (
                ev_per_mvgd.drop(columns=["bus_id"])
                .groupby(["scenario", "scenario_variation"])
                .count()
            )

            for (
                scn,
                scn_var,
            ), ev_actual in ev_per_mvgd_counts_per_scn.iterrows():
                scenario_parameters = get_sector_parameters(
                    "mobility", scenario=scn
                )["motorized_individual_travel"]

                # Get EV target
                ev_target = scenario_parameters[scn_var]["ev_count"]

                np.testing.assert_allclose(
                    int(ev_actual),
                    ev_target,
                    rtol=0.0001,
                    err_msg=f"Dataset on EV numbers allocated to MVGDs "
                    f"seems to be flawed. "
                    f"Scenario: [{scn}], "
                    f"Scenario variation: [{scn_var}].",
                )


1			"""
2			* Calculate number of electric vehicles and allocate on different spatial
3			levels: :py:func:`allocate_evs_numbers`
4			* Allocate specific EVs to MV grid districts:
5			:py:func:`allocate_evs_to_grid_districts`
6
7			"""
8
9			from itertools import permutations
10
11			from sqlalchemy.sql import func
12			import numpy as np
13			import pandas as pd
14
15			from egon.data import db
16			from egon.data.datasets.emobility.motorized_individual_travel.db_classes import (
17			EgonEvCountMunicipality,
18			EgonEvCountMvGridDistrict,
19			EgonEvCountRegistrationDistrict,
20			EgonEvMvGridDistrict,
21			EgonEvPool,
22			)
23			from egon.data.datasets.emobility.motorized_individual_travel.helpers import (
24			COLUMNS_KBA,
25			CONFIG_EV,
26			TESTMODE_OFF,
27			read_kba_data,
28			read_rs7_data,
29			)
30			from egon.data.datasets.emobility.motorized_individual_travel.tests import (
31			validate_electric_vehicles_numbers,
32			)
33			from egon.data.datasets.mv_grid_districts import MvGridDistricts
34			from egon.data.datasets.scenario_parameters import get_sector_parameters
35			from egon.data.datasets.zensus_mv_grid_districts import MapZensusGridDistricts
36			from egon.data.datasets.zensus_vg250 import (
37			DestatisZensusPopulationPerHaInsideGermany,
38			MapZensusVg250,
39			Vg250Gem,
40			Vg250GemPopulation,
41			)
42			import egon.data.config
43
44			RANDOM_SEED = egon.data.config.settings()["egon-data"]["--random-seed"]
45
46
47			def fix_missing_ags_municipality_regiostar(muns, rs7_data):
48			"""Check if all AGS of municipality dataset are included in RegioStaR7
49			dataset and vice versa.
50
51			As of Dec 2021, some municipalities are not included int the RegioStaR7
52			dataset. This is mostly caused by incorporations of a municipality by
53			another municipality. This is fixed by assigning a RS7 id from another
54			municipality with similar AGS (most likely a neighboured one).
55
56			Missing entries in the municipality dataset is printed but not fixed
57			as it doesn't result in bad data. Nevertheless, consider to update the
58			municipality/VG250 dataset.
59
60			Parameters
61			----------
62			muns : pandas.DataFrame
63			Municipality data
64			rs7_data : pandas.DataFrame
65			RegioStaR7 data
66
67			Returns
68			-------
69			pandas.DataFrame
70			Fixed RegioStaR7 data
71			"""
72
73			if len(muns.ags.unique()) != len(rs7_data.ags):
74			print(
75			"==========> Number of AGS differ between VG250 and RS7, "
76			"trying to fix this..."
77			)
78
79			# Get AGS differences
80			ags_datasets = {"RS7": rs7_data.ags, "VG250": muns.ags}
81			ags_datasets_missing = {k: [] for k in ags_datasets.keys()}
82			perm = permutations(ags_datasets.items())
83			for (name1, ags1), (name2, ags2) in perm:
84			print(f" Checking if all AGS of {name1} are in {name2}...")
85			missing = [_ for _ in ags1 if _ not in ags2.to_list()]
86			if len(missing) > 0:
87			ags_datasets_missing[name2] = missing
88			print(f" AGS in {name1} but not in {name2}: ", missing)
89			else:
90			print(" OK")
91
92			print("")
93
94			# Try to fix
95			for name, missing in ags_datasets_missing.items():
96			if len(missing) > 0:
97			# RS7 entries missing: use RS7 number from mun with similar AGS
98			if name == "RS7":
99			for ags in missing:
100			similar_entry = rs7_data[
101			round((rs7_data.ags.div(10))) == round(ags / 10)
102			].iloc[0]
103			if len(similar_entry) > 0:
104			print(
105			f"Adding dataset from VG250 to RS7 "
106			f"based upon AGS {ags}."
107			)
108			similar_entry.ags = ags
109			rs7_data = rs7_data.append(similar_entry)
110			print("Consider to update RS7.")
111			# VG250 entries missing:
112			elif name == "VG250":
113			print(
114			"Cannot guess VG250 entries. This error does not "
115			"result in bad data but consider to update VG250."
116			)
117
118			if len(muns.ags.unique()) != len(rs7_data.ags):
119			print("==========> AGS could not be fixed!")
120			else:
121			print("==========> AGS were fixed!")
122
123			return rs7_data
124
125
126			def calc_evs_per_reg_district(scenario_variation_parameters, kba_data):
127			"""Calculate EVs per registration district
128
129			Parameters
130			----------
131			scenario_variation_parameters : dict
132			Parameters of scenario variation
133			kba_data : pandas.DataFrame
134			Vehicle registration data for registration district
135
136			Returns
137			-------
138			pandas.DataFrame
139			EVs per registration district
140			"""
141
142			scenario_variation_parameters["mini_share"] = (
143			scenario_variation_parameters["bev_mini_share"]
144			+ scenario_variation_parameters["phev_mini_share"]
145			)
146			scenario_variation_parameters["medium_share"] = (
147			scenario_variation_parameters["bev_medium_share"]
148			+ scenario_variation_parameters["phev_medium_share"]
149			)
150			scenario_variation_parameters["luxury_share"] = (
151			scenario_variation_parameters["bev_luxury_share"]
152			+ scenario_variation_parameters["phev_luxury_share"]
153			)
154
155			factor_dict = dict()
156			factor_dict["mini_factor"] = (
157			scenario_variation_parameters["mini_share"]
158			* scenario_variation_parameters["ev_count"]
159			/ kba_data.mini.sum()
160			)
161			factor_dict["medium_factor"] = (
162			scenario_variation_parameters["medium_share"]
163			* scenario_variation_parameters["ev_count"]
164			/ kba_data.medium.sum()
165			)
166			factor_dict["luxury_factor"] = (
167			scenario_variation_parameters["luxury_share"]
168			* scenario_variation_parameters["ev_count"]
169			/ kba_data.luxury.sum()
170			)
171
172			# Define shares and factors
173			ev_data = kba_data.copy()
174
175			for tech, params in CONFIG_EV.items():
176			ev_data[tech] = (
177			(
178			kba_data[params["column"]]
179			* factor_dict[params["factor"]]
180			* scenario_variation_parameters[params["tech_share"]]
181			/ scenario_variation_parameters[params["share"]]
182			)
183			.round()
184			.astype("int")
185			)
186
187			ev_data.drop(
188			columns=[_ for _ in COLUMNS_KBA if _ != "reg_district"], inplace=True
189			)
190
191			return ev_data
192
193
194			def calc_evs_per_municipality(ev_data, rs7_data):
195			"""Calculate EVs per municipality
196
197			Parameters
198			----------
199			ev_data : pandas.DataFrame
200			EVs per regstration district
201			rs7_data : pandas.DataFrame
202			RegioStaR7 data
203			"""
204			with db.session_scope() as session:
205			query = session.query(
206			Vg250GemPopulation.ags_0.label("ags"),
207			Vg250GemPopulation.gen,
208			Vg250GemPopulation.population_total.label("pop"),
209			)
210
211			muns = pd.read_sql(
212			query.statement, query.session.bind, index_col=None
213			).astype({"ags": "int64"})
214
215			muns["ags_district"] = (
216			muns.ags.multiply(1 / 1000).apply(np.floor).astype("int")
217			)
218
219			# Manual fix of Trier-Saarburg: Introduce new `ags_reg_district`
220			# for correct allocation of mun to registration district
221			# (Zulassungsbezirk), see above for background.
222			muns["ags_reg_district"] = muns["ags_district"]
223			muns.loc[muns["ags_reg_district"] == 7235, "ags_reg_district"] = 7211
224
225			# Remove multiple municipality entries (due to 'gf' in VG250)
226			# by summing up population
227			muns = (
228			muns[["ags", "gen", "ags_reg_district", "pop"]]
229			.groupby(["ags", "gen", "ags_reg_district"])
230			.sum()
231			.reset_index()
232			)
233
234			# Add population of registration district
235			pop_per_reg_district = (
236			muns[["ags_reg_district", "pop"]]
237			.groupby("ags_reg_district")
238			.sum()
239			.rename(columns={"pop": "pop_district"})
240			.reset_index()
241			)
242
243			# Fix missing ags in mun data if not in testmode
244			if TESTMODE_OFF:
245			rs7_data = fix_missing_ags_municipality_regiostar(muns, rs7_data)
246
247			# Merge municipality, EV data and pop per district
248			ev_data_muns = muns.merge(ev_data, on="ags_reg_district").merge(
249			pop_per_reg_district, on="ags_reg_district"
250			)
251
252			# Disaggregate EV numbers to municipality
253			for tech in ev_data[CONFIG_EV.keys()]:
254			ev_data_muns[tech] = round(
255			ev_data_muns[tech]
256			* ev_data_muns["pop"]
257			/ ev_data_muns["pop_district"]
258			).astype("int")
259
260			# Filter columns
261			cols = ["ags"]
262			cols.extend(CONFIG_EV.keys())
263			ev_data_muns = ev_data_muns[cols]
264
265			# Merge RS7 data
266			ev_data_muns = ev_data_muns.merge(rs7_data[["ags", "rs7_id"]], on="ags")
267
268			return ev_data_muns
269
270
271			def calc_evs_per_grid_district(ev_data_muns):
272			"""Calculate EVs per grid district by using population weighting
273
274			Parameters
275			----------
276			ev_data_muns : pandas.DataFrame
277			EV data for municipalities
278
279			Returns
280			-------
281			pandas.DataFrame
282			EV data for grid districts
283			"""
284
285			# Read MVGDs with intersecting muns and aggregate pop for each
286			# municipality part
287			with db.session_scope() as session:
288			query_pop_per_mvgd = (
289			session.query(
290			MvGridDistricts.bus_id,
291			Vg250Gem.ags,
292			func.sum(
293			DestatisZensusPopulationPerHaInsideGermany.population
294			).label("pop"),
295			)
296			.select_from(MapZensusGridDistricts)
297			.join(
298			MvGridDistricts,
299			MapZensusGridDistricts.bus_id == MvGridDistricts.bus_id,
300			)
301			.join(
302			DestatisZensusPopulationPerHaInsideGermany,
303			MapZensusGridDistricts.zensus_population_id
304			== DestatisZensusPopulationPerHaInsideGermany.id,
305			)
306			.join(
307			MapZensusVg250,
308			MapZensusGridDistricts.zensus_population_id
309			== MapZensusVg250.zensus_population_id,
310			)
311			.join(
312			Vg250Gem, MapZensusVg250.vg250_municipality_id == Vg250Gem.id
313			)
314			.group_by(MvGridDistricts.bus_id, Vg250Gem.ags)
315			.order_by(Vg250Gem.ags)
316			)
317
318			mvgd_pop_per_mun = pd.read_sql(
319			query_pop_per_mvgd.statement,
320			query_pop_per_mvgd.session.bind,
321			index_col=None,
322			).astype({"bus_id": "int64", "pop": "int64", "ags": "int64"})
323
324			# Calc population share of each municipality in MVGD
325			mvgd_pop_per_mun_in_mvgd = mvgd_pop_per_mun.groupby(["bus_id", "ags"]).agg(
326			{"pop": "sum"}
327			)
328
329			# Calc relative and absolute population shares:
330			# * pop_mun_in_mvgd: pop share of mun which intersects with MVGD
331			# * pop_share_mun_in_mvgd: relative pop share of mun which
332			# intersects with MVGD
333			# * pop_mun_total: total pop of mun
334			# * pop_mun_in_mvgd_of_mun_total: relative pop share of mun which
335			# intersects with MVGD in relation to total pop of mun
336			mvgd_pop_per_mun_in_mvgd = (
337			mvgd_pop_per_mun_in_mvgd.groupby(level=0)
338			.apply(lambda x: x / float(x.sum()))
339			.reset_index()
340			.rename(columns={"pop": "pop_share_mun_in_mvgd"})
341			.merge(
342			mvgd_pop_per_mun_in_mvgd.reset_index(),
343			on=["bus_id", "ags"],
344			how="left",
345			)
346			.rename(columns={"pop": "pop_mun_in_mvgd"})
347			.merge(
348			mvgd_pop_per_mun[["ags", "pop"]]
349			.groupby("ags")
350			.agg({"pop": "sum"}),
351			on="ags",
352			how="left",
353			)
354			.rename(columns={"pop": "pop_mun_total"})
355			)
356			mvgd_pop_per_mun_in_mvgd["pop_mun_in_mvgd_of_mun_total"] = (
357			mvgd_pop_per_mun_in_mvgd["pop_mun_in_mvgd"]
358			/ mvgd_pop_per_mun_in_mvgd["pop_mun_total"]
359			)
360
361			# Merge EV data
362			ev_data_mvgds = mvgd_pop_per_mun_in_mvgd.merge(
363			ev_data_muns, on="ags", how="left"
364			).sort_values(["bus_id", "ags"])
365
366			# Calc EVs per MVGD by using EV from mun and share of mun's pop
367			# that is located within MVGD
368			for tech in ev_data_mvgds[CONFIG_EV.keys()]:
369			ev_data_mvgds[tech] = (
370			round(
371			ev_data_mvgds[tech]
372			* ev_data_mvgds["pop_mun_in_mvgd_of_mun_total"]
373			)
374			.fillna(0)
375			.astype("int")
376			)
377
378			# Set RS7 id for MVGD by using the RS7 id from the mun with the
379			# highest share in population
380			rs7_data_mvgds = (
381			ev_data_mvgds[["bus_id", "pop_mun_in_mvgd", "rs7_id"]]
382			.groupby(["bus_id", "rs7_id"])
383			.sum()
384			.sort_values(
385			["bus_id", "pop_mun_in_mvgd"], ascending=False, na_position="last"
386			)
387			.reset_index()
388			.drop_duplicates("bus_id", keep="first")[["bus_id", "rs7_id"]]
389			)
390
391			# Join RS7 id and select columns
392			columns = ["bus_id"] + [_ for _ in CONFIG_EV.keys()]
393			ev_data_mvgds = (
394			ev_data_mvgds[columns]
395			.groupby("bus_id")
396			.agg("sum")
397			.merge(rs7_data_mvgds, on="bus_id", how="left")
398			)
399
400			return ev_data_mvgds
401
402
403			def allocate_evs_numbers():
404			"""Allocate electric vehicles to different spatial levels.
405
406			Accocation uses today's vehicles registration data per registration
407			district from KBA and scales scenario's EV targets (BEV and PHEV)
408			linearly using population. Furthermore, a RegioStaR7 code (BMVI) is
409			assigned.
410
411			Levels:
412			* districts of registration
413			* municipalities
414			* grid districts
415
416			Parameters
417			----------
418
419			Returns
420			-------
421
422			"""
423			# Import
424			kba_data = read_kba_data()
425			rs7_data = read_rs7_data()
426
427			for scenario_name in ["eGon2035", "eGon100RE"]:
428			# Load scenario params
429			scenario_parameters = get_sector_parameters(
430			"mobility", scenario=scenario_name
431			)["motorized_individual_travel"]
432
433			print(f"SCENARIO: {scenario_name}")
434
435			# Go through scenario variations
436			for (
437			scenario_variation_name,
438			scenario_variation_parameters,
439			) in scenario_parameters.items():
440
441			print(f" SCENARIO VARIATION: {scenario_variation_name}")
442
443			# Get EV target
444			ev_target = scenario_variation_parameters["ev_count"]
445
446			#####################################
447			# EV data per registration district #
448			#####################################
449			print("Calculate EV numbers for registration districts...")
450			ev_data = calc_evs_per_reg_district(
451			scenario_variation_parameters, kba_data
452			)
453			# Check EV results if not in testmode
454			if TESTMODE_OFF:
455			validate_electric_vehicles_numbers(
456			"EV count in registration districts", ev_data, ev_target
457			)
458			# Add scenario columns and write to DB
459			ev_data["scenario"] = scenario_name
460			ev_data["scenario_variation"] = scenario_variation_name
461			ev_data.sort_values(
462			["scenario", "scenario_variation", "ags_reg_district"],
463			inplace=True,
464			)
465			ev_data.to_sql(
466			name=EgonEvCountRegistrationDistrict.__table__.name,
467			schema=EgonEvCountRegistrationDistrict.__table__.schema,
468			con=db.engine(),
469			if_exists="append",
470			index=False,
471			)
472
473			#####################################
474			# EV data per municipality #
475			#####################################
476			print("Calculate EV numbers for municipalities...")
477			ev_data_muns = calc_evs_per_municipality(ev_data, rs7_data)
478			# Check EV results if not in testmode
479			if TESTMODE_OFF:
480			validate_electric_vehicles_numbers(
481			"EV count in municipalities", ev_data_muns, ev_target
482			)
483			# Add scenario columns and write to DB
484			ev_data_muns["scenario"] = scenario_name
485			ev_data_muns["scenario_variation"] = scenario_variation_name
486			ev_data_muns.sort_values(
487			["scenario", "scenario_variation", "ags"], inplace=True
488			)
489			ev_data_muns.to_sql(
490			name=EgonEvCountMunicipality.__table__.name,
491			schema=EgonEvCountMunicipality.__table__.schema,
492			con=db.engine(),
493			if_exists="append",
494			index=False,
495			)
496
497			#####################################
498			# EV data per grid district #
499			#####################################
500			print("Calculate EV numbers for grid districts...")
501			ev_data_mvgds = calc_evs_per_grid_district(ev_data_muns)
502			# Check EV results if not in testmode
503			if TESTMODE_OFF:
504			validate_electric_vehicles_numbers(
505			"EV count in grid districts", ev_data_mvgds, ev_target
506			)
507			# Add scenario columns and write to DB
508			ev_data_mvgds["scenario"] = scenario_name
509			ev_data_mvgds["scenario_variation"] = scenario_variation_name
510			ev_data_mvgds.sort_values(
511			["scenario", "scenario_variation", "bus_id"], inplace=True
512			)
513			ev_data_mvgds.to_sql(
514			name=EgonEvCountMvGridDistrict.__table__.name,
515			schema=EgonEvCountMvGridDistrict.__table__.schema,
516			con=db.engine(),
517			if_exists="append",
518			index=False,
519			)
520
521
522			def allocate_evs_to_grid_districts():
523			"""Allocate EVs to MV grid districts for all scenarios and scenario
524			variations.
525
526			Each grid district in
527			:class:`egon.data.datasets.mv_grid_districts.MvGridDistricts`
528			is assigned a list of electric vehicles from the EV pool in
529			:class:`EgonEvPool` based on the RegioStar7 region and the
530			counts per EV type in :class:`EgonEvCountMvGridDistrict`.
531			Results are written to :class:`EgonEvMvGridDistrict`.
532			"""
533
534			def get_random_evs(row):
535			"""Get random EV sample for EV type and RS7 region"""
536			return (
537			ev_pool.loc[
			0 ignored issues – show introduced 2022-06-03 12:34 UTC by Report Bug Copy Issue Report The variable `ev_pool` does not seem to be defined in case the `for` loop on line `545` is not entered. Are you sure this can never be the case? Loading history...
538			(ev_pool.rs7_id == row.rs7_id)
539			& (ev_pool["type"] == row["type"])
540			]
541			.sample(row["count"], replace=True)
542			.ev_id.to_list()
543			)
544
545			for scenario_name in ["eGon2035", "eGon100RE"]:
546			print(f"SCENARIO: {scenario_name}")
547
548			# Load EVs per grid district
549			print("Loading EV counts for grid districts...")
550			with db.session_scope() as session:
551			query = session.query(EgonEvCountMvGridDistrict).filter(
552			EgonEvCountMvGridDistrict.scenario == scenario_name
553			)
554			ev_per_mvgd = pd.read_sql(
555			query.statement, query.session.bind, index_col=None
556			)
557
558			# Convert EV types' wide to long format
559			ev_per_mvgd = pd.melt(
560			ev_per_mvgd,
561			id_vars=["scenario", "scenario_variation", "bus_id", "rs7_id"],
562			value_vars=CONFIG_EV.keys(),
563			var_name="type",
564			value_name="count",
565			)
566
567			# Load EV pool
568			print(" Loading EV pool...")
569			with db.session_scope() as session:
570			query = session.query(EgonEvPool).filter(
571			EgonEvPool.scenario == scenario_name
572			)
573			ev_pool = pd.read_sql(
574			query.statement,
575			query.session.bind,
576			index_col=None,
577			)
578
579			# Draw EVs randomly for each grid district from pool
580			print(" Draw EVs from pool for grid districts...")
581			np.random.seed(RANDOM_SEED)
582			ev_per_mvgd["egon_ev_pool_ev_id"] = ev_per_mvgd.apply(
583			get_random_evs, axis=1
584			)
585			ev_per_mvgd.drop(columns=["rs7_id", "type", "count"], inplace=True)
586
587			# EV lists to rows
588			ev_per_mvgd = ev_per_mvgd.explode("egon_ev_pool_ev_id")
589
590			# Check for empty entries
591			empty_ev_entries = ev_per_mvgd.egon_ev_pool_ev_id.isna().sum()
592			if empty_ev_entries > 0:
593			print("====================================================")
594			print(
595			f"WARNING: Found {empty_ev_entries} empty entries "
596			f"and will remove it:"
597			)
598			print(ev_per_mvgd[ev_per_mvgd.egon_ev_pool_ev_id.isna()])
599			ev_per_mvgd = ev_per_mvgd[~ev_per_mvgd.egon_ev_pool_ev_id.isna()]
600			print("====================================================")
601
602			# Write trips to DB
603			print(" Writing allocated data to DB...")
604			ev_per_mvgd.to_sql(
605			name=EgonEvMvGridDistrict.__table__.name,
606			schema=EgonEvMvGridDistrict.__table__.schema,
607			con=db.engine(),
608			if_exists="append",
609			index=False,
610			method="multi",
611			chunksize=10000,
612			)
613
614			# Check EV result sums for all scenario variations if not in testmode
615			if TESTMODE_OFF:
616			print(" Validating results...")
617			ev_per_mvgd_counts_per_scn = (
618			ev_per_mvgd.drop(columns=["bus_id"])
619			.groupby(["scenario", "scenario_variation"])
620			.count()
621			)
622
623			for (
624			scn,
625			scn_var,
626			), ev_actual in ev_per_mvgd_counts_per_scn.iterrows():
627			scenario_parameters = get_sector_parameters(
628			"mobility", scenario=scn
629			)["motorized_individual_travel"]
630
631			# Get EV target
632			ev_target = scenario_parameters[scn_var]["ev_count"]
633
634			np.testing.assert_allclose(
635			int(ev_actual),
636			ev_target,
637			rtol=0.0001,
638			err_msg=f"Dataset on EV numbers allocated to MVGDs "
639			f"seems to be flawed. "
640			f"Scenario: [{scn}], "
641			f"Scenario variation: [{scn_var}].",
642			)
643

openego / eGon-data

Push — dev ( 85d654...312d71 )

calc_evs_per_grid_district() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like