Merge pull request #367 from openego/bugfix/cts_timeseries

Bugfix/cts timeseries

edisgo/flex_opt/charging_strategies.py

@@ -128,7 +128,7 @@ def charging_strategy(
             f"to the original frequency of the edisgo time series data."
         )
 
-        edisgo_obj.resample_timeseries(freq=simbev_timedelta)
+        edisgo_obj.timeseries.resample_timeseries(freq=simbev_timedelta)
 
     if strategy == "dumb":
         # "dumb" charging
@@ -309,7 +309,7 @@ def charging_strategy(
         raise ValueError(f"Strategy {strategy} has not yet been implemented.")
 
     if resample:
-        edisgo_obj.resample_timeseries(freq=edisgo_timedelta)
+        edisgo_obj.timeseries.resample_timeseries(freq=edisgo_timedelta)
 
     # set reactive power time series to 0 Mvar
     # fmt: off

edisgo/io/dsm_import.py

@@ -215,6 +215,16 @@ def get_profile_cts(
         hourly resolution in MW. Index contains hour of the year (from 0 to 8759) and
         column names are site ID as integer.
 
+    Notes
+    ------
+    Be aware, that in this function the DSM time series are disaggregated to all CTS
+    loads in the grid. In some cases, this can lead to an over- or underestimation of
+    the DSM potential, as in egon_data buildings are mapped to a grid based on the
+    zensus cell they are in whereas in ding0 buildings are mapped to a grid based on
+    the geolocation. As it can happen that buildings lie outside an MV grid but within
+    a zensus cell that is assigned to that MV grid, they are mapped differently in
+    egon_data and ding0.
+
     """
     saio.register_schema("demand", engine)
     from saio.demand import egon_etrago_electricity_cts_dsm_timeseries

edisgo/io/generators_import.py

@@ -817,6 +817,15 @@ def oedb(
     engine : :sqlalchemy:`sqlalchemy.Engine<sqlalchemy.engine.Engine>`
         Database engine.
 
+    Notes
+    ------
+    Note, that PV rooftop plants are queried using the building IDs not the MV grid ID
+    as in egon_data buildings are mapped to a grid based on the
+    zensus cell they are in whereas in ding0 buildings are mapped to a grid based on
+    the geolocation. As it can happen that buildings lie outside an MV grid but within
+    a zensus cell that is assigned to that MV grid, they are mapped differently in
+    egon_data and ding0, and it is therefore better to query using the building IDs.
+
     """
 
     def _get_egon_power_plants():
@@ -853,14 +862,15 @@ def _get_egon_power_plants():
             subtype=power_plants_gdf["type"].map(mapping)
         )
         # unwrap source ID
-        power_plants_gdf["source_id"] = power_plants_gdf.apply(
-            lambda _: (
-                list(_["source_id"].values())[0]
-                if isinstance(_["source_id"], dict)
-                else None
-            ),
-            axis=1,
-        )
+        if not power_plants_gdf.empty:
+            power_plants_gdf["source_id"] = power_plants_gdf.apply(
+                lambda _: (
+                    list(_["source_id"].values())[0]
+                    if isinstance(_["source_id"], dict)
+                    else None
+                ),
+                axis=1,
+            )
         return power_plants_gdf
 
     def _get_egon_pv_rooftop():
@@ -935,6 +945,33 @@ def _get_egon_chp_plants():
     power_plants_gdf = _get_egon_power_plants()
     chp_gdf = _get_egon_chp_plants()
 
+    # sanity check - kick out generators that are too large
+    p_nom_upper = edisgo_object.config["grid_connection"]["upper_limit_voltage_level_4"]
+    drop_ind = pv_rooftop_df[pv_rooftop_df.p_nom > p_nom_upper].index
+    if len(drop_ind) > 0:
+        logger.warning(
+            f"There are {len(drop_ind)} PV rooftop plants with a nominal capacity "
+            f"larger {p_nom_upper} MW. Connecting them to the MV is not valid, "
+            f"wherefore they are dropped."
+        )
+        pv_rooftop_df.drop(index=drop_ind, inplace=True)
+    drop_ind = power_plants_gdf[power_plants_gdf.p_nom > p_nom_upper].index
+    if len(drop_ind) > 0:
+        logger.warning(
+            f"There are {len(drop_ind)} power plants with a nominal capacity "
+            f"larger {p_nom_upper} MW. Connecting them to the MV is not valid, "
+            f"wherefore they are dropped."
+        )
+        power_plants_gdf.drop(index=drop_ind, inplace=True)
+    drop_ind = chp_gdf[chp_gdf.p_nom > p_nom_upper].index
+    if len(drop_ind) > 0:
+        logger.warning(
+            f"There are {len(drop_ind)} CHP plants with a nominal capacity "
+            f"larger {p_nom_upper} MW. Connecting them to the MV is not valid, "
+            f"wherefore they are dropped."
+        )
+        chp_gdf.drop(index=drop_ind, inplace=True)
+
     # determine number of generators and installed capacity in future scenario
     # for validation of grid integration
     total_p_nom_scenario = (

edisgo/io/heat_pump_import.py

@@ -126,8 +126,13 @@ def _get_central_heat_pumps():
                 ),
             )
         )
-
-        df = gpd.read_postgis(query.statement, engine, index_col=None)
+        srid = db.get_srid_of_db_table(session, egon_district_heating.geometry)
+        df = gpd.read_postgis(
+            query.statement,
+            engine,
+            index_col=None,
+            crs=f"EPSG:{srid}",
+        )
 
         # transform to same SRID as MV grid district geometry
         return df.to_crs(mv_grid_geom_srid)

edisgo/io/timeseries_import.py

@@ -437,7 +437,7 @@ def heat_demand_oedb(edisgo_obj, scenario, engine, timeindex=None):
             building_ids, scenario, engine
         )
         cts_profiles_df = get_cts_profiles_per_building(
-            edisgo_obj.topology.id, scenario, "heat", engine
+            edisgo_obj, scenario, "heat", engine
         )
         # drop CTS profiles for buildings without a heat pump
         buildings_no_hp = [_ for _ in cts_profiles_df.columns if _ not in building_ids]
@@ -462,7 +462,10 @@ def heat_demand_oedb(edisgo_obj, scenario, engine, timeindex=None):
         individual_heating_df = pd.DataFrame(index=timeindex_full)
 
     # get district heating profiles from oedb
-    dh_ids = hp_df.district_heating_id.dropna().unique()
+    if "district_heating_id" in hp_df.columns:
+        dh_ids = hp_df.district_heating_id.dropna().unique()
+    else:
+        dh_ids = []
     if len(dh_ids) > 0:
         dh_profile_df = get_district_heating_heat_demand_profiles(
             dh_ids, scenario, engine
@@ -575,7 +578,7 @@ def electricity_demand_oedb(
     cts_building_ids = cts_loads.building_id.dropna().unique()
     if len(cts_building_ids) > 0:
         cts_profiles_df = get_cts_profiles_per_building(
-            edisgo_obj.topology.id, scenario, "electricity", engine
+            edisgo_obj, scenario, "electricity", engine
         )
         drop_buildings = [
             _ for _ in cts_profiles_df.columns if _ not in cts_building_ids
@@ -932,14 +935,78 @@ def get_district_heating_heat_demand_profiles(district_heating_ids, scenario, en
     return df.astype("float")
 
 
-def get_cts_profiles_per_building(
+def get_cts_profiles_per_building(edisgo_obj, scenario, sector, engine):
+    """
+    Gets CTS heat demand profiles per CTS building for all CTS buildings in MV grid.
+
+    This function is a helper function that should not be but is necessary, as in
+    egon_data buildings are mapped to a grid based on the zensus cell they are in
+    whereas in ding0 buildings are mapped to a grid based on the geolocation. As it can
+    happen that buildings lie outside an MV grid but within a zensus cell that is
+    assigned to that MV grid, they are mapped differently in egon_data and ding0.
+    This function therefore checks, if there are CTS loads with other grid IDs and if
+    so, gets profiles for other grid IDs (by calling
+    :func:`~.io.timeseries_import.get_cts_profiles_per_grid` with different grid IDs)
+    in order to obtain a demand profile for all CTS loads.
+
+    Parameters
+    ----------
+    edisgo_obj : :class:`~.EDisGo`
+    scenario : str
+        Scenario for which to retrieve demand data. Possible options
+        are 'eGon2035' and 'eGon100RE'.
+    sector : str
+        Demand sector for which profile is calculated: "electricity" or "heat"
+    engine : :sqlalchemy:`sqlalchemy.Engine<sqlalchemy.engine.Engine>`
+        Database engine.
+
+    Returns
+    -------
+    :pandas:`pandas.DataFrame<DataFrame>`
+        Dataframe with CTS demand profiles per building for one year in an
+        hourly resolution in MW. Index contains hour of the year (from 0 to 8759) and
+        column names are building ID as integer.
+
+    """
+    saio.register_schema("boundaries", engine)
+    from saio.boundaries import egon_map_zensus_mvgd_buildings
+
+    # get MV grid IDs of CTS loads
+    cts_loads = edisgo_obj.topology.loads_df[
+        (edisgo_obj.topology.loads_df.type == "conventional_load")
+        & (edisgo_obj.topology.loads_df.sector == "cts")
+    ]
+    cts_building_ids = cts_loads.building_id.dropna().unique()
+    with session_scope_egon_data(engine) as session:
+        query = session.query(
+            egon_map_zensus_mvgd_buildings.building_id,
+            egon_map_zensus_mvgd_buildings.bus_id,
+        ).filter(
+            egon_map_zensus_mvgd_buildings.building_id.in_(cts_building_ids),
+        )
+        df = pd.read_sql(query.statement, engine, index_col="building_id")
+
+    # iterate over grid IDs
+    profiles_df = pd.DataFrame()
+    for bus_id in df.bus_id.unique():
+        profiles_grid_df = get_cts_profiles_per_grid(
+            bus_id=bus_id, scenario=scenario, sector=sector, engine=engine
+        )
+        profiles_df = pd.concat([profiles_df, profiles_grid_df], axis=1)
+
+    # filter CTS loads in grid
+    return profiles_df.loc[:, cts_building_ids]
+
+
+def get_cts_profiles_per_grid(
     bus_id,
     scenario,
     sector,
     engine,
 ):
     """
-    Gets CTS heat demand profiles per building for all buildings in MV grid.
+    Gets CTS heat or electricity demand profiles per building for all buildings in the
+    given MV grid.
 
     Parameters
     ----------
@@ -951,7 +1018,7 @@ def get_cts_profiles_per_building(
     sector : str
         Demand sector for which profile is calculated: "electricity" or "heat"
     engine : :sqlalchemy:`sqlalchemy.Engine<sqlalchemy.engine.Engine>`
-            Database engine.
+        Database engine.
 
     Returns
     -------
@@ -1252,7 +1319,7 @@ def _get_profiles(profile_ids):
 
     # calculate demand profile per building
     ts_df = pd.DataFrame()
-    for building_id, df in profile_ids_buildings.groupby(by=["building_id"]):
+    for building_id, df in profile_ids_buildings.groupby(by="building_id"):
         load_ts_building = (
             profiles_df.loc[:, df["profile_id"]].sum(axis=1)
             * df["factor"].iloc[0]

edisgo/network/timeseries.py

@@ -2190,7 +2190,7 @@ def resample_timeseries(
                 self.timeindex[0],
                 self.timeindex[-1] + freq_orig,
                 freq=freq,
-                closed="left",
+                inclusive="left",
             )
         else:  # down-sampling
             index = pd.date_range(

rtd_requirements.txt

@@ -8,7 +8,7 @@ jupyter_dash
 matplotlib >= 3.3.0
 multiprocess
 networkx >= 2.5.0
-pandas >= 1.2.0
+pandas >= 1.2.0, < 2.0.0
 plotly
 pyomo >= 6.0
 pypower

setup.py

@@ -45,7 +45,7 @@ def read(fname):
     "matplotlib >= 3.3.0",
     "multiprocess",
     "networkx >= 2.5.0",
-    "pandas >= 1.2.0",
+    "pandas >= 1.2.0, < 2.0.0",
     "plotly",
     "pydot",
     "pygeos",

tests/io/test_pypsa_io.py

@@ -2,7 +2,7 @@
 import pandas as pd
 import pytest
 
-from pandas.util.testing import assert_frame_equal
+from pandas.testing import assert_frame_equal
 
 from edisgo import EDisGo
 from edisgo.io import pypsa_io

tests/io/test_timeseries_import.py

@@ -228,11 +228,33 @@ def test_get_district_heating_heat_demand_profiles(self):
     @pytest.mark.local
     def test_get_cts_profiles_per_building(self):
 
+        edisgo_object = EDisGo(
+            ding0_grid=pytest.ding0_test_network_3_path, legacy_ding0_grids=False
+        )
+        cts_loads = edisgo_object.topology.loads_df[
+            edisgo_object.topology.loads_df.sector == "cts"
+        ]
+        df = timeseries_import.get_cts_profiles_per_building(
+            edisgo_object, "eGon2035", "electricity", pytest.engine
+        )
+        assert df.shape == (8760, len(cts_loads))
+
+        # manipulate CTS load to lie within another grid
+        edisgo_object.topology.loads_df.at[cts_loads.index[0], "building_id"] = 5
         df = timeseries_import.get_cts_profiles_per_building(
+            edisgo_object, "eGon2035", "electricity", pytest.engine
+        )
+        assert df.shape == (8760, len(cts_loads))
+        # ToDo add further tests
+
+    @pytest.mark.local
+    def test_get_cts_profiles_per_grid(self):
+
+        df = timeseries_import.get_cts_profiles_per_grid(
             33535, "eGon2035", "heat", pytest.engine
         )
         assert df.shape == (8760, 85)
-        df = timeseries_import.get_cts_profiles_per_building(
+        df = timeseries_import.get_cts_profiles_per_grid(
             33535, "eGon2035", "electricity", pytest.engine
         )
         assert df.shape == (8760, 85)

tests/network/test_electromobility.py

@@ -8,7 +8,7 @@
 import pandas as pd
 import pytest
 
-from pandas.util.testing import assert_frame_equal
+from pandas.testing import assert_frame_equal
 
 from edisgo.edisgo import EDisGo
 from edisgo.io import electromobility_import

tests/network/test_timeseries.py

@@ -8,11 +8,7 @@
 import pandas as pd
 import pytest
 
-from pandas.util.testing import (
-    assert_frame_equal,
-    assert_index_equal,
-    assert_series_equal,
-)
+from pandas.testing import assert_frame_equal, assert_index_equal, assert_series_equal
 
 from edisgo import EDisGo
 from edisgo.network import timeseries

tests/test_edisgo.py

@@ -11,7 +11,7 @@
 import pytest
 
 from matplotlib import pyplot as plt
-from pandas.util.testing import assert_frame_equal, assert_series_equal
+from pandas.testing import assert_frame_equal, assert_series_equal
 from shapely.geometry import Point
 
 from edisgo import EDisGo