Hotfix for #193 (#194)

edisgo now loads all weather cells that intersect with mv grid

Co-authored-by: birgits <birgit.schachler@rl-institut.de>

edisgo/network/timeseries.py

@@ -6,10 +6,10 @@
 
 from workalendar.europe import Germany
 from demandlib import bdew as bdew, particular_profiles as profiles
-
 from edisgo.io.timeseries_import import import_feedin_timeseries
 from edisgo.tools.tools import assign_voltage_level_to_component,\
-    drop_duplicated_columns
+    drop_duplicated_columns, get_weather_cells_intersecting_with_grid_district
+
 
 logger = logging.getLogger("edisgo")
 
@@ -597,7 +597,6 @@ def get_component_timeseries(edisgo_obj, **kwargs):
         ranges of the given time series that will be used in the analysis.
 
     """
-
     mode = kwargs.get("mode", None)
     timeindex = kwargs.get("timeindex", edisgo_obj.timeseries.timeindex)
     # reset TimeSeries
@@ -647,17 +646,19 @@ def get_component_timeseries(edisgo_obj, **kwargs):
             raise ValueError("{} is not a valid mode.".format(mode))
     else:
         config_data = edisgo_obj.config
-        weather_cell_ids = (
-            edisgo_obj.topology.generators_df.weather_cell_id.dropna().unique()
-        )
+
+        weather_cell_ids = get_weather_cells_intersecting_with_grid_district(
+            edisgo_obj)
+
         # feed-in time series of fluctuating renewables
         ts = kwargs.get("timeseries_generation_fluctuating", None)
         if isinstance(ts, pd.DataFrame):
             edisgo_obj.timeseries.generation_fluctuating = ts
         elif isinstance(ts, str) and ts == "oedb":
-            edisgo_obj.timeseries.generation_fluctuating = import_feedin_timeseries(
-                config_data, weather_cell_ids, kwargs.get("timeindex", None)
-            )
+            edisgo_obj.timeseries.generation_fluctuating = \
+                import_feedin_timeseries(
+                    config_data, weather_cell_ids, kwargs.get(
+                        "timeindex", None))
         else:
             raise ValueError(
                 "Your input for "

edisgo/tools/tools.py

@@ -1,10 +1,28 @@
+import os
+import numpy as np
 import pandas as pd
 import networkx as nx
 from math import pi, sqrt
+from sqlalchemy import func
 
 from edisgo.flex_opt import exceptions
 from edisgo.flex_opt import check_tech_constraints
 from edisgo.network.grids import LVGrid
+from edisgo.tools import session_scope
+
+if "READTHEDOCS" not in os.environ:
+
+    from egoio.db_tables import climate
+    from egoio.tools.db import connection
+
+    from shapely.geometry.multipolygon import MultiPolygon
+    from shapely.wkt import loads as wkt_loads
+
+    geopandas = True
+    try:
+        import geopandas as gpd
+    except:
+        geopandas = False
 
 
 def select_worstcase_snapshots(edisgo_obj):
@@ -406,3 +424,58 @@ def assign_voltage_level_to_component(edisgo_obj, df):
         axis=1,
     )
     return df
+
+
+def get_weather_cells_intersecting_with_grid_district(edisgo_obj):
+    """
+    Get all weather cells that intersect with the grid district.
+    
+    Parameters
+    ----------
+    edisgo_obj : :class:`~.EDisGo`
+
+    Returns
+    -------
+    set
+        Set with weather cell IDs
+
+    """
+
+    # Download geometries of weather cells
+    srid = edisgo_obj.topology.grid_district["srid"]
+    table = climate.Cosmoclmgrid
+    with session_scope() as session:
+        query = session.query(
+            table.gid,
+            func.ST_AsText(
+                func.ST_Transform(
+                    table.geom, srid
+                )
+            ).label("geometry")
+        )
+    geom_data = pd.read_sql_query(
+        query.statement, query.session.bind)
+    geom_data.geometry = geom_data.apply(
+        lambda _: wkt_loads(_.geometry), axis=1)
+    geom_data = gpd.GeoDataFrame(
+        geom_data, crs=f"EPSG:{srid}")
+
+    # Make sure MV Geometry is MultiPolygon
+    mv_geom = edisgo_obj.topology.grid_district["geom"]
+    if mv_geom.geom_type == "Polygon":
+        # Transform Polygon to MultiPolygon and overwrite geometry
+        p = wkt_loads(str(mv_geom))
+        m = MultiPolygon([p])
+        edisgo_obj.topology.grid_district["geom"] = m
+    elif mv_geom.geom_type == "MultiPolygon":
+        m = mv_geom
+    else:
+        raise ValueError(
+            f"Grid district geometry is of type {type(mv_geom)}."
+            " Only Shapely Polygon or MultiPolygon are accepted.")
+    mv_geom_gdf = gpd.GeoDataFrame(
+        m, crs=f"EPSG:{srid}", columns=["geometry"])
+
+    return set(np.append(gpd.sjoin(
+        geom_data, mv_geom_gdf, how="right", op='intersects').gid.unique(),
+        edisgo_obj.topology.generators_df.weather_cell_id.dropna().unique()))

setup.py

@@ -58,7 +58,9 @@ def create_edisgo_path():
         'egoio >= 0.4.7',
         'matplotlib >= 3.3.0',
         'pypower',
-        'sklearn'
+        'sklearn',
+        'pydot',
+        'Rtree',
     ],
     extras_require={
         'geoplot': ['geopandas >= 0.9.0', 'contextily', 'descartes'],

tests/tools/test_tools.py

@@ -191,3 +191,16 @@ def test_assign_feeder(self):
         assert not topo.buses_df[
             ~topo.buses_df.index.isin(
                 mv_buses)].lv_feeder.isna().any()
+
+    def test_get_weather_cells_intersecting_with_grid_district(self):
+        weather_cells = \
+            tools.get_weather_cells_intersecting_with_grid_district(
+                self.edisgo)
+        assert len(weather_cells) == 4
+        assert 1123075 in weather_cells
+        assert 1122075 in weather_cells
+        assert 1122076 in weather_cells
+        # the following weather cell does not intersect with the grid district
+        # but there are generators in the grid that have that weather cell
+        # for some reason..
+        assert 1122074 in weather_cells