diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e6e4b7b..5a8552b4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ All notable changes to this project will be documented in this file. - Specific `Line(s)Result` and line utilisation calculation [#176](https://github.com/ie3-institute/pypsdm/issues/176) - Big refactoring of result types to extract more generic time series [#192](https://github.com/ie3-institute/pypsdm/pull/192) - Result types do not contain uuid and optional name anymore, but time series dicts now have Entity key that contain the information [#192](https://github.com/ie3-institute/pypsdm/pull/192) +- Add WeatherDict data types and retrieval of weighted nodal weather [#193](https://github.com/ie3-institute/pypsdm/issues/193) ### Changed diff --git a/pypsdm/db/weather/models.py b/pypsdm/db/weather/models.py index 2642313d..d3866f49 100644 --- a/pypsdm/db/weather/models.py +++ b/pypsdm/db/weather/models.py @@ -47,6 +47,7 @@ def wind_velocity_v(self): def name_mapping(): return { "time": "time", + "coordinate_id": "coordinate_id", "aswdifd_s": "diffuse_irradiance", "aswdir_s": "direct_irradiance", "t2m": "temperature", @@ -56,7 +57,7 @@ def name_mapping(): class Coordinate(SQLModel, table=True): - id: Optional[int] = Field(default=None, primary_key=True) + id: int = Field(primary_key=True) coordinate: str = Field() def __eq__(self, other): @@ -86,3 +87,8 @@ def longitude(self) -> float: @property def x(self) -> float: return self.point.x + + @staticmethod + def from_xy(id, x, y): + wkb = Point(x, y).wkb_hex + return Coordinate(id=id, coordinate=wkb) diff --git a/pypsdm/db/weather/proxy.py b/pypsdm/db/weather/proxy.py index 7a9648e3..53e13338 100644 --- a/pypsdm/db/weather/proxy.py +++ b/pypsdm/db/weather/proxy.py @@ -68,7 +68,7 @@ def get_closest_coordinates( self, x: float, y: float, - n: int, + n: int, # amount of closest coordinates to return schema_name="public", table_name="coordinate", id_column="id", @@ -125,74 +125,3 @@ def create_engine_from_params( return create_engine( f"postgresql://{username}:{password}@{host}:{port}/{database}", echo=echo ) - - -def weighted_interpolation_coordinates( - target: tuple[float, float], - nearest_coords: list[tuple[Coordinate, float]], -) -> list[tuple[Coordinate, float]]: - """ - Given a list of nearest surrounding cordinates with respect to a target coordinate, - find the nearest coordinate in each quadrant and weigh them by their distance to - the target. - - Requires at least one coordinate in each quadrant (meaing top left, top right, - bottom left, bottom right). - - Args: - target (tuple[float, float]): Target coordinate (x (longitude), y (latitude)) - nearest_coords (list[tuple[Coordinate, float]]): List of nearest coordinates - with their distances to the target - """ - - x, y = target - - # Check if the queried coordinate is surrounded in each quadrant - quadrants: list[tuple[Coordinate | None, float]] = [ - (None, float("inf")) for _ in range(4) - ] # [Q1, Q2, Q3, Q4] - for point, distance in nearest_coords: - - if point.x < x and point.y > y: - if quadrants[0][0]: - if distance < quadrants[0][1]: - quadrants[0] = (point, distance) - else: - quadrants[0] = (point, distance) - - if point.x > x and point.y > y: - if quadrants[1][0]: - if distance < quadrants[1][1]: - quadrants[1] = (point, distance) - else: - quadrants[1] = (point, distance) - - if point.x < x and point.y < y: - if quadrants[2][0]: - if distance < quadrants[2][1]: - quadrants[2] = (point, distance) - else: - quadrants[2] = (point, distance) - - if point.x > x and point.y < y: - if quadrants[3][0]: - if distance < quadrants[3][1]: - quadrants[3] = (point, distance) - else: - quadrants[3] = (point, distance) - - acc_dist = 0 - for q in quadrants: - if q[0]: - acc_dist += q[1] - else: - raise ValueError("Not all quadrants are filled") - - n = len(quadrants) - weighted_coordinates = [] - for q in quadrants: - if q[0]: - weight = (1 - (q[1] / acc_dist)) / (n - 1) - weighted_coordinates.append((q[0], weight)) - - return weighted_coordinates diff --git a/pypsdm/db/weather/utils.py b/pypsdm/db/weather/utils.py new file mode 100644 index 00000000..29108671 --- /dev/null +++ b/pypsdm/db/weather/utils.py @@ -0,0 +1,141 @@ +import concurrent.futures +from datetime import datetime + +from pypsdm.db.weather.models import Coordinate +from pypsdm.db.weather.proxy import WeatherProxy +from pypsdm.models.input.node import Nodes +from pypsdm.models.ts.types import CoordinateWeather, WeatherDict + + +def get_nodal_weighted_weather( + nodes: Nodes, start: datetime, end: datetime, weather: WeatherProxy +) -> WeatherDict: + weighted_coordinates = nodal_weighted_coordinates(nodes=nodes, weather=weather) + coordinates = set() + for wc in weighted_coordinates.values(): + coordinates.update([c[0].id for c in wc]) + values = weather.get_weather_for_interval(start, end, coordinates) + weather_dct = WeatherDict.from_value_list(values) + + weighted_weather_dct = {} + for node_id, wc in weighted_coordinates.items(): + weighted_weather = CoordinateWeather.empty() + for c, w in wc: + weighted_weather += weather_dct[c.id] * w + weighted_weather_dct[node_id] = weighted_weather + + return WeatherDict(weighted_weather_dct) + + +def nodal_weighted_coordinates( + nodes: Nodes, weather: WeatherProxy +) -> dict[str, list[tuple[Coordinate, float]]]: + """ + Determine the 4 nearest coordinates (each in one of the surrounding quadrants) for + each of the nodes and weigh them by their distance (sum of the weights = 1). + + Args: + nodes (Nodes): Nodes object containing the nodes + weather (WeatherProxy): WeatherProxy object to access the weather database + + Returns: + dict[str, list[tuple[Coordinate, float]]]: Dictionary with node uuids as keys + and a list of tuples with the nearest coordinates and their weights as values + """ + nodes_uuid = nodes.uuid + weighted_coordinates = {} + + def fetch_and_weight_coordinates(node: str, nodes: Nodes, weather: WeatherProxy): + lon, lat = nodes.data.loc[node, ["longitude", "latitude"]] # type: ignore + closest = weather.get_closest_coordinates(lon, lat, 8) + weighted_coord = weighted_interpolation_coordinates((lon, lat), closest) + return node, weighted_coord + + with concurrent.futures.ThreadPoolExecutor(max_workers=10) as executor: + future_to_node = { + executor.submit(fetch_and_weight_coordinates, node, nodes, weather): node + for node in nodes_uuid + } + + for future in concurrent.futures.as_completed(future_to_node): + node = future_to_node[future] + try: + node, weighted_coord = future.result() + weighted_coordinates[node] = weighted_coord + except Exception as exc: + print(f"Node {node} generated an exception: {exc}") + return weighted_coordinates + + +def weighted_interpolation_coordinates( + target: tuple[float, float], + nearest_coords: list[tuple[Coordinate, float]], +) -> list[tuple[Coordinate, float]]: + """ + Given a list of nearest surrounding cordinates with respect to a target coordinate, + find the nearest coordinate in each quadrant and weigh them by their distance to + the target. + + Requires at least one coordinate in each quadrant (meaing top left, top right, + bottom left, bottom right). + + Note that the nearest coordinates can be found with the find n closest functionality + of the WeatherProxy. + + Args: + target (tuple[float, float]): Target coordinate (x (longitude), y (latitude)) + nearest_coords (list[tuple[Coordinate, float]]): List of nearest coordinates + with their distances to the target + """ + + x, y = target + + # Check if the queried coordinate is surrounded in each quadrant + quadrants: list[tuple[Coordinate | None, float]] = [ + (None, float("inf")) for _ in range(4) + ] # [Q1, Q2, Q3, Q4] + for point, distance in nearest_coords: + + if point.x < x and point.y > y: + if quadrants[0][0]: + if distance < quadrants[0][1]: + quadrants[0] = (point, distance) + else: + quadrants[0] = (point, distance) + + if point.x > x and point.y > y: + if quadrants[1][0]: + if distance < quadrants[1][1]: + quadrants[1] = (point, distance) + else: + quadrants[1] = (point, distance) + + if point.x < x and point.y < y: + if quadrants[2][0]: + if distance < quadrants[2][1]: + quadrants[2] = (point, distance) + else: + quadrants[2] = (point, distance) + + if point.x > x and point.y < y: + if quadrants[3][0]: + if distance < quadrants[3][1]: + quadrants[3] = (point, distance) + else: + quadrants[3] = (point, distance) + + acc_dist = 0 + for q in quadrants: + if q[0]: + acc_dist += q[1] + else: + raise ValueError("Not all quadrants are filled") + + n = len(quadrants) + weighted_coordinates = [] + for q in quadrants: + if q[0]: + weight = (1 - (q[1] / acc_dist)) / (n - 1) + weighted_coordinates.append((q[0], weight)) + + return weighted_coordinates diff --git a/pypsdm/models/result/participant/dict.py b/pypsdm/models/result/participant/dict.py index fa1d98e0..dc6bd9d4 100644 --- a/pypsdm/models/result/participant/dict.py +++ b/pypsdm/models/result/participant/dict.py @@ -25,6 +25,8 @@ class EntitiesResultDictMixin: + def uuids(self) -> set[str]: + return {key.uuid for key in self.keys()} # type: ignore @classmethod @abstractmethod @@ -142,7 +144,6 @@ def from_csv_for_entity( class EmsResult(ComplexPowerDict[EntityKey], EntitiesResultDictMixin): - def __init__(self, data: dict[EntityKey, ComplexPower]): for key, value in data.items(): if not isinstance(key, EntityKey): @@ -167,7 +168,6 @@ def entity_type(cls) -> EntitiesEnum: class LoadsResult(ComplexPowerDict[EntityKey], EntitiesResultDictMixin): - def __init__(self, data: dict[EntityKey, ComplexPower]): for key, value in data.items(): if not isinstance(key, EntityKey): @@ -192,7 +192,6 @@ def entity_type(cls) -> EntitiesEnum: class FixedFeedInsResult(ComplexPowerDict[EntityKey], EntitiesResultDictMixin): - def __init__(self, data: dict[EntityKey, ComplexPower]): for key, value in data.items(): if not isinstance(key, EntityKey): @@ -217,7 +216,6 @@ def entity_type(cls) -> EntitiesEnum: class PvsResult(ComplexPowerDict[EntityKey], EntitiesResultDictMixin): - def __init__(self, data: dict[EntityKey, ComplexPower]): for key, value in data.items(): if not isinstance(key, EntityKey): @@ -242,7 +240,6 @@ def entity_type(cls) -> EntitiesEnum: class WecsResult(ComplexPowerDict[EntityKey], EntitiesResultDictMixin): - def __init__(self, data: dict[EntityKey, ComplexPower]): for key, value in data.items(): if not isinstance(key, EntityKey): @@ -267,7 +264,6 @@ def entity_type(cls) -> EntitiesEnum: class StoragesResult(ComplexPowerWithSocDict[EntityKey], EntitiesResultDictMixin): - def __init__(self, data: dict[EntityKey, ComplexPowerWithSoc]): for key, value in data.items(): if not isinstance(key, EntityKey): diff --git a/pypsdm/models/ts/mixins.py b/pypsdm/models/ts/mixins.py index b5c8d989..d5a43985 100644 --- a/pypsdm/models/ts/mixins.py +++ b/pypsdm/models/ts/mixins.py @@ -210,3 +210,44 @@ def v_complex( ) -> DataFrame: # ffill not supported for complex numbers return self.attr_df("v_complex", False, favor_ids, v_rated) + + +class WeatherDataMixin(AttributeMixin): + + @property + def diffuse_irradiance(self): + return self.data["diffuse_irradiance"] # type: ignore + + @property + def direct_irradiance(self): + return self.data["direct_irradiance"] # type: ignore + + @property + def temperature(self): + return self.data["temperature"] # type: ignore + + @property + def wind_velocity_u(self): + return self.data["wind_velocity_u"] # type: ignore + + @property + def wind_velocity_v(self): + return self.data["wind_velocity_v"] # type: ignore + + +class WeatherDataDictMixin(TimeSeriesDictMixin): + + def diffuse_irradiance(self, ffill=True, favor_ids=True) -> DataFrame: + return self.attr_df("diffuse_irradiance", ffill, favor_ids) + + def direct_irradiance(self, ffill=True, favor_ids=True) -> DataFrame: + return self.attr_df("direct_irradiance", ffill, favor_ids) + + def temperature(self, ffill=True, favor_ids=True) -> DataFrame: + return self.attr_df("temperature", ffill, favor_ids) + + def wind_velocity_u(self, ffill=True, favor_ids=True) -> DataFrame: + return self.attr_df("wind_velocity_u", ffill, favor_ids) + + def wind_velocity_v(self, ffill=True, favor_ids=True) -> DataFrame: + return self.attr_df("wind_velocity_v", ffill, favor_ids) diff --git a/pypsdm/models/ts/types.py b/pypsdm/models/ts/types.py index ef92d163..af5b0b55 100644 --- a/pypsdm/models/ts/types.py +++ b/pypsdm/models/ts/types.py @@ -1,5 +1,9 @@ +from __future__ import annotations + from functools import reduce -from typing import List, Self, Sequence, Union +from typing import TYPE_CHECKING, List, Self, Sequence, Union + +from pandas import DataFrame from pypsdm.models.ts.base import K, TimeSeries, TimeSeriesDict from pypsdm.models.ts.mixins import ( @@ -10,10 +14,15 @@ ComplexVoltageMixin, SocDictMixin, SocMixin, + WeatherDataDictMixin, + WeatherDataMixin, ) from pypsdm.processing.dataframe import add_df from pypsdm.processing.series import Tuple, add_series +if TYPE_CHECKING: + from pypsdm.db.weather.models import WeatherValue + class ComplexPower(TimeSeries, ComplexPowerMixin): def __eq__(self, other: object) -> bool: @@ -170,3 +179,55 @@ def __eq__(self, other: object) -> bool: def complex_power_sum(self) -> ComplexPower: return ComplexPower.sum(list(self.values())) + + +class CoordinateWeather(TimeSeries, WeatherDataMixin): + + def __add__(self, other) -> Self: + return self.__class__(add_df(self.data, other.data)) + + def __mul__(self, other: float | int) -> Self: + return self.__class__(self.data * other) + + __rmul__ = __mul__ + + @classmethod + def from_value_list(cls, values: list[WeatherValue]): + df = cls.df_from_value_list(values) + if df["coordinate_id"].nunique() > 1: + raise ValueError("Multiple coordinate ids in weather data") + df.drop(columns=["coordinate_id"], inplace=True) + return cls(df) + + @staticmethod + def df_from_value_list(values: list[WeatherValue]): + from pypsdm.db.weather.models import WeatherValue + + value_dicts = [value.model_dump() for value in values] + df = DataFrame(value_dicts, columns=WeatherValue.__table__.columns.keys()) + df.rename(columns=WeatherValue.name_mapping(), inplace=True) + df.set_index("time", inplace=True, drop=True) + return df + + @staticmethod + def attributes() -> List[str]: + return [ + "diffuse_irradiance", + "direct_irradiance", + "temperature", + "wind_velocity_u", + "wind_velocity_v", + ] + + +class WeatherDict(TimeSeriesDict[int, CoordinateWeather], WeatherDataDictMixin): + + @classmethod + def from_value_list(cls, values: list[WeatherValue]): + df = CoordinateWeather.df_from_value_list(values) + grps = df.groupby("coordinate_id") + dct = {} + for coord_id, grp in grps: + grp = grp.drop(columns=["coordinate_id"]) + dct[coord_id] = CoordinateWeather(grp) + return cls(dct) diff --git a/pypsdm/processing/dataframe.py b/pypsdm/processing/dataframe.py index 48087e68..edf8cfca 100644 --- a/pypsdm/processing/dataframe.py +++ b/pypsdm/processing/dataframe.py @@ -11,7 +11,10 @@ def add_df(a: pd.DataFrame, b: pd.DataFrame): Adds two dataframes with different indices in an event discrete manner. """ if not a.columns.equals(b.columns): - raise ValueError("DataFrames have different columns") + diff = set(a.columns).symmetric_difference(set(b.columns)) + raise ValueError( + f"DataFrames have different columns: {diff} not in both DataFrames." + ) if len(a) == 0: return b diff --git a/tests/db/weather/test_proxy.py b/tests/db/weather/test_proxy.py index b128caf8..31acdce6 100644 --- a/tests/db/weather/test_proxy.py +++ b/tests/db/weather/test_proxy.py @@ -1,7 +1,7 @@ import math from pypsdm.db.weather.models import Coordinate -from pypsdm.db.weather.proxy import weighted_interpolation_coordinates +from pypsdm.db.weather.utils import weighted_interpolation_coordinates def test_weighted_interpolation_coordinates(): diff --git a/tests/db/weather/test_utils.py b/tests/db/weather/test_utils.py new file mode 100644 index 00000000..e874245a --- /dev/null +++ b/tests/db/weather/test_utils.py @@ -0,0 +1,36 @@ +from pypsdm.db.weather.models import Coordinate +from pypsdm.db.weather.utils import weighted_interpolation_coordinates + + +def test_weighted_interpolation_coordinates(): + target = (0, 0) + + coord1 = Coordinate.from_xy(id=1, x=1, y=1) + coord2 = Coordinate.from_xy(id=2, x=-1, y=1) + coord3 = Coordinate.from_xy(id=3, x=1, y=-1) + coord4 = Coordinate.from_xy(id=4, x=-1, y=-1) + coord5 = Coordinate.from_xy(id=5, x=-2, y=-2) + + # Calculate distances manually for the test setup + nearest_coords = [ + (coord2, 1.414), + (coord1, 1.414), + (coord4, 1.414), + (coord3, 1.414), + (coord5, 2.828), + ] + + result = weighted_interpolation_coordinates(target, nearest_coords) + + expected_weight = 0.25 + expected = [ + (coord2, expected_weight), + (coord1, expected_weight), + (coord4, expected_weight), + (coord3, expected_weight), + ] + + # Verify that the function returns the expected results + assert sorted(result, key=lambda x: x[0].id) == sorted( + expected, key=lambda x: x[0].id + ), "Weights or coordinates do not match expected values" diff --git a/tests/models/ts/test_types.py b/tests/models/ts/test_types.py index e73e6350..a70d0c81 100644 --- a/tests/models/ts/test_types.py +++ b/tests/models/ts/test_types.py @@ -1,5 +1,8 @@ +from datetime import datetime + import pandas as pd +from pypsdm.db.weather.models import WeatherValue from pypsdm.models.ts.base import TIME_COLUMN_NAME, EntityKey from pypsdm.models.ts.types import ( ComplexPower, @@ -8,6 +11,8 @@ ComplexVoltage, ComplexVoltageDict, ComplexVoltagePower, + CoordinateWeather, + WeatherDict, ) @@ -178,3 +183,55 @@ def test_voltage_power(): vp = get_voltage_power_data() v = vp.as_complex_voltage() pd.testing.assert_series_equal(vp.v_mag, v.v_mag) + + +def test_weather_from_value_list(): + weather_value1 = WeatherValue( + time=datetime(2021, 1, 1, 0, 0), + coordinate_id=1, + aswdifd_s=1, + aswdir_s=2, + t2m=3, + u131m=4, + v131m=5, + ) + weather_value2 = WeatherValue( + time=datetime(2021, 1, 1, 2, 0), + coordinate_id=1, + aswdifd_s=6, + aswdir_s=7, + t2m=8, + u131m=9, + v131m=10, + ) + + values = [weather_value1, weather_value2] + df = CoordinateWeather.df_from_value_list(values) + expected_cols = set(CoordinateWeather.attributes()) + expected_cols.add("coordinate_id") + assert set(df.columns) == expected_cols + assert set(df.index) == set([weather_value1.time, weather_value2.time]) + + weather_value2.coordinate_id = 2 + dct = WeatherDict.from_value_list(values) + assert set(dct.keys()) == set([1, 2]) + + +def test_weather_add_mult(): + data = pd.DataFrame( + { + "diffuse_irradiance": [1, 2, 3], + "direct_irradiance": [4, 5, 6], + "temperature": [7, 8, 9], + "wind_velocity_u": [10, 11, 12], + "wind_velocity_v": [13, 14, 15], + }, + index=pd.date_range("2021-01-01", periods=3, freq="H"), + ) + data.index.name = TIME_COLUMN_NAME + weather = CoordinateWeather(data) + weather2 = CoordinateWeather(data) + res = weather + weather2 + pd.testing.assert_frame_equal(res.data, data * 2, check_dtype=False) + res = weather + CoordinateWeather.empty() + pd.testing.assert_frame_equal(res.data, data, check_dtype=False)