Create waking_state.py

wind_up/waking_state.py

@@ -0,0 +1,294 @@
+import numpy as np
+import pandas as pd
+from geographiclib.geodesic import Geodesic
+from tabulate import tabulate
+
+from wind_up.constants import (
+    DEFAULT_AIR_DENSITY,
+    RAW_DOWNTIME_S_COL,
+    RAW_POWER_COL,
+    RAW_WINDSPEED_COL,
+    TIMEBASE_S,
+    TIMESTAMP_COL,
+)
+from wind_up.math_funcs import circ_diff
+from wind_up.models import PlotConfig, TurbineType, WindUpConfig
+from wind_up.plots.waking_state_plots import plot_waking_state_one_ttype_or_wtg
+from wind_up.wind_funcs import calc_cp
+
+
+def add_waking_state_one_ttype(wf_df: pd.DataFrame, ttype: TurbineType, plot_cfg: PlotConfig | None) -> pd.DataFrame:
+    wf_df = wf_df.copy()
+    rated_power = ttype.rated_power_kw
+    wf_df["waking"] = ~wf_df["ActivePowerMean"].isna()
+    margin_for_waking_rated_power = 0.8
+    wf_df["waking"] = wf_df["waking"] | (wf_df[RAW_POWER_COL] > margin_for_waking_rated_power * rated_power)
+    margin_for_cp_calc = 0.3
+    df_for_median_cp = wf_df[
+        wf_df["waking"]
+        & (wf_df[RAW_POWER_COL] > margin_for_cp_calc * rated_power)
+        & (wf_df[RAW_POWER_COL] < (1 - margin_for_cp_calc) * rated_power)
+    ]
+    median_waking_cp = calc_cp(  # type: ignore[union-attr]
+        df_for_median_cp["ActivePowerMean"],
+        df_for_median_cp["WindSpeedMean"].clip(lower=1),
+        DEFAULT_AIR_DENSITY,
+        ttype.rotor_diameter_m,
+    ).median()
+    factor_for_waking_cp = 0.7
+    wf_df["waking"] = wf_df["waking"] | (
+        calc_cp(
+            wf_df[RAW_POWER_COL],
+            wf_df[RAW_WINDSPEED_COL].clip(lower=1),
+            DEFAULT_AIR_DENSITY,
+            ttype.rotor_diameter_m,
+        )
+        > factor_for_waking_cp * median_waking_cp
+    )
+
+    margin_for_not_waking_rated_power = 0.01
+    wf_df["not_waking"] = wf_df[RAW_POWER_COL] < margin_for_not_waking_rated_power * rated_power
+    wf_df["not_waking"] = wf_df["not_waking"] | (wf_df[RAW_DOWNTIME_S_COL] > TIMEBASE_S * factor_for_waking_cp)
+
+    wf_df["not_waking"] = wf_df["not_waking"] & (~wf_df["waking"])
+    wf_df["unknown_waking"] = (~wf_df["waking"]) & (~wf_df["not_waking"])
+
+    if plot_cfg is not None:
+        waking_frc = wf_df["waking"].sum() / len(wf_df)
+        print(f"{ttype.turbine_type} {waking_frc * 100:.1f}% of rows are waking")
+        not_waking_frc = wf_df["not_waking"].sum() / len(wf_df)
+        print(
+            f"{ttype.turbine_type} {not_waking_frc * 100:.1f}% of rows are not waking",
+        )
+        unknown_waking_frc = wf_df["unknown_waking"].sum() / len(wf_df)
+        print(
+            f"{ttype.turbine_type} {unknown_waking_frc * 100:.1f}% of rows have unknown or partial waking",
+        )
+        plot_waking_state_one_ttype_or_wtg(wf_df=wf_df, ttype_or_wtg=ttype.turbine_type, plot_cfg=plot_cfg)
+        for wtg_name in wf_df.index.unique(level="TurbineName"):
+            df_wtg = wf_df.loc[wtg_name]
+            plot_waking_state_one_ttype_or_wtg(wf_df=df_wtg, ttype_or_wtg=wtg_name, plot_cfg=plot_cfg)
+    return wf_df
+
+
+def add_waking_state(cfg: WindUpConfig, wf_df: pd.DataFrame, plot_cfg: PlotConfig | None) -> pd.DataFrame:
+    df_input = wf_df.copy()
+    wf_df = pd.DataFrame()
+    for ttype in cfg.list_unique_turbine_types():
+        wtgs = cfg.list_turbine_ids_of_type(ttype)
+        df_ttype = df_input.loc[wtgs]
+        wf_df_ = add_waking_state_one_ttype(
+            wf_df=df_ttype,
+            ttype=ttype,
+            plot_cfg=plot_cfg,
+        )
+        wf_df = pd.concat([wf_df, wf_df_])
+    return wf_df.sort_index()
+
+
+def calc_bearing(*, lat1: float, long1: float, lat2: float, long2: float) -> float:
+    bearing_deg = Geodesic.WGS84.Inverse(lat1, long1, lat2, long2)["azi1"]
+    return bearing_deg % 360
+
+
+def calc_distance(*, lat1: float, long1: float, lat2: float, long2: float) -> float:
+    return Geodesic.WGS84.Inverse(lat1, long1, lat2, long2)["s12"]
+
+
+distance_and_bearing_cache: dict[tuple[float, float, float, float], tuple[float, float]] = {}
+
+
+def get_distance_and_bearing(*, lat1: float, long1: float, lat2: float, long2: float) -> tuple[float, float]:
+    # see if distance and bearing are in known_answers
+    # if so, return them
+    if (lat1, long1, lat2, long2) in distance_and_bearing_cache:
+        distance_m, bearing_deg = distance_and_bearing_cache[(lat1, long1, lat2, long2)]
+    else:
+        distance_m = calc_distance(lat1=lat1, long1=long1, lat2=lat2, long2=long2)
+        bearing_deg = calc_bearing(lat1=lat1, long1=long1, lat2=lat2, long2=long2)
+        distance_and_bearing_cache[(lat1, long1, lat2, long2)] = (distance_m, bearing_deg)
+    return distance_m, bearing_deg
+
+
+def calc_iec_upwind_turbines(*, lat: float, long: float, wind_direction: float, cfg: WindUpConfig) -> list[str]:
+    # find all turbines in cfg which are upwind of lat, long
+    wind_direction = wind_direction % 360
+    wtg_names = []
+    bearings = []
+    distances = []
+    wtg_d_norms = []
+    for wtg in cfg.asset.wtgs:
+        wtg_lat = wtg.latitude
+        wtg_long = wtg.longitude
+        wtg_names.append(wtg.name)
+        wtg_distance, wtg_bearing = get_distance_and_bearing(lat1=lat, long1=long, lat2=wtg_lat, long2=wtg_long)
+        bearings.append(wtg_bearing)
+        distances.append(wtg_distance)
+        wtg_d_norms.append(wtg_distance / wtg.turbine_type.rotor_diameter_m)
+    upwind_df = pd.DataFrame(
+        data={
+            "wtg_name": wtg_names,
+            "bearing": bearings,
+            "distance_m": distances,
+            "distance_diameters": wtg_d_norms,
+        },
+    )
+    upwind_df["disturbed_sector"] = 180.0
+    ge_2_diameters = upwind_df["distance_diameters"] >= 2  # noqa PLR2004
+    upwind_df.loc[ge_2_diameters, "disturbed_sector"] = (
+        1.3 * np.rad2deg(np.arctan(2.5 / upwind_df.loc[ge_2_diameters, "distance_diameters"] + 0.15)) + 10
+    )
+    gt_20_diameters = upwind_df["distance_diameters"] > 20  # noqa PLR2004
+    upwind_df.loc[gt_20_diameters, "disturbed_sector"] = 0
+
+    upwind_df["relative_bearing"] = circ_diff(upwind_df["bearing"], wind_direction)
+    upwind_df["iec_upwind"] = abs(upwind_df["relative_bearing"]) < upwind_df["disturbed_sector"] / 2
+
+    upwind_turbine_list = list(upwind_df.query("iec_upwind")["wtg_name"].sort_values().values)
+    if wind_direction % 90 == 0:
+        print(
+            f"calc_iec_upwind_turbines lat={lat:.2f} long={long:.2f} "
+            f"wind_dir={wind_direction:.0f} {upwind_turbine_list}",
+        )
+    return upwind_turbine_list
+
+
+upwind_wtgs_cache: dict[tuple[float, float, float, str | None], list[str]] = {}
+
+
+def lat_long_is_valid(lat: float, long: float) -> bool:
+    return (-90 <= lat <= 90) and (-180 <= long <= 180)  # noqa PLR2004
+
+
+def get_iec_upwind_turbines_one_latlong(
+    *,
+    lat: float,
+    long: float,
+    wind_direction: float,
+    cfg: WindUpConfig,
+    object_name: str | None = None,
+) -> list[str]:
+    if not lat_long_is_valid(lat, long):
+        msg = f"lat={lat} long={long} is not a valid lat long"
+        raise ValueError(msg)
+
+    if (lat, long, wind_direction, object_name) in upwind_wtgs_cache:
+        upwind_wtgs = upwind_wtgs_cache[(lat, long, wind_direction, object_name)]
+    else:
+        upwind_wtgs = calc_iec_upwind_turbines(lat=lat, long=long, wind_direction=wind_direction, cfg=cfg)
+        if object_name is not None:
+            upwind_wtgs = [x for x in upwind_wtgs if x.lower() != object_name.lower()]
+        upwind_wtgs_cache[(lat, long, wind_direction, object_name)] = upwind_wtgs
+    return upwind_wtgs
+
+
+def get_iec_upwind_turbines(
+    *,
+    latlongs: list[tuple[float, float]],
+    wind_direction: float,
+    cfg: WindUpConfig,
+    object_name: str | None = None,
+) -> list[str]:
+    upwind_wtgs = []
+    for lat, long in latlongs:
+        upwind_wtgs += get_iec_upwind_turbines_one_latlong(
+            lat=lat,
+            long=long,
+            wind_direction=wind_direction,
+            cfg=cfg,
+            object_name=object_name,
+        )
+    return sorted(set(upwind_wtgs))
+
+
+def calc_scen_name_from_wtgs_not_waking(wtgs_not_waking: list[str]) -> str:
+    if len(wtgs_not_waking) < 1:
+        msg = "wtgs_not_waking must have at least one element"
+        raise ValueError(msg)
+    return " ".join(wtgs_not_waking) + " offline"
+
+
+def list_wtgs_offline_in_scen(waking_scenario: str) -> list[str]:
+    return [x for x in waking_scenario.split(" ") if x != "offline"]
+
+
+def add_waking_scen(
+    *,
+    test_name: str,
+    ref_name: str,
+    test_ref_df: pd.DataFrame,
+    cfg: WindUpConfig,
+    wf_df: pd.DataFrame,
+    ref_wd_col: str,
+    ref_lat: float,
+    ref_long: float,
+) -> pd.DataFrame:
+    test_ref_df = test_ref_df.copy()
+    test_ref_df["waking_scenario"] = "not calculated"
+    try:
+        test_wtg = next(x for x in cfg.asset.wtgs if x.name == test_name)
+    except StopIteration as exc:
+        msg = f"{test_name} not found in cfg.asset.wtgs"
+        raise ValueError(msg) from exc
+
+    all_turbines_waking = wf_df.groupby(TIMESTAMP_COL)["waking"].all().to_frame(name="all_turbines_waking")
+    test_ref_df = test_ref_df.merge(all_turbines_waking, how="left", left_index=True, right_index=True)
+    test_ref_df.loc[test_ref_df["all_turbines_waking"], "waking_scenario"] = "none offline"
+
+    test_ref_df["rounded_wd"] = test_ref_df[ref_wd_col].round(0).mod(360)
+    waking_df = wf_df[["waking", "not_waking", "unknown_waking"]].pivot_table(
+        index=TIMESTAMP_COL,
+        columns="TurbineName",
+        aggfunc="max",
+        observed=False,
+    )
+
+    for wd in sorted(test_ref_df["rounded_wd"].unique()):
+        if wd >= 0:
+            test_upwind_wtgs = get_iec_upwind_turbines(
+                latlongs=test_wtg.get_latlongs(),
+                wind_direction=wd,
+                cfg=cfg,
+                object_name=test_name,
+            )
+            if lat_long_is_valid(ref_lat, ref_long):
+                ref_upwind_wtgs = get_iec_upwind_turbines(
+                    latlongs=[(ref_lat, ref_long)],
+                    wind_direction=wd,
+                    cfg=cfg,
+                    object_name=ref_name,
+                )
+            else:
+                ref_upwind_wtgs = []
+            upwind_wtgs = sorted(set(test_upwind_wtgs + ref_upwind_wtgs))
+            if len(upwind_wtgs) == 0:
+                test_ref_df.loc[test_ref_df["rounded_wd"] == wd, "waking_scenario"] = "none offline"
+            else:
+                for idx in test_ref_df[
+                    (test_ref_df["rounded_wd"] == wd) & (test_ref_df["waking_scenario"] == "not calculated")
+                ].index:
+                    if waking_df.loc[idx, "unknown_waking"].loc[upwind_wtgs].any():
+                        test_ref_df.loc[idx, "waking_scenario"] = "unknown"
+                    elif waking_df.loc[idx, "waking"].loc[upwind_wtgs].all():
+                        test_ref_df.loc[idx, "waking_scenario"] = "none offline"
+                    else:
+                        not_waking_series = waking_df.loc[idx, "not_waking"].loc[upwind_wtgs]
+                        test_ref_df.loc[idx, "waking_scenario"] = calc_scen_name_from_wtgs_not_waking(
+                            list(not_waking_series[not_waking_series].index),
+                        )
+
+    if (test_ref_df.dropna(subset=ref_wd_col)["waking_scenario"] == "not calculated").any():
+        msg = "some rows with defined ref_wd_col have waking scenario = 'not calculated'"
+        raise RuntimeError(msg)
+
+    top_scens = test_ref_df.dropna(subset=ref_wd_col)["waking_scenario"].value_counts()[0:5]
+    print(f"top {len(top_scens)} {test_name} {ref_name} waking scenarios [%]:")
+    print(
+        tabulate(
+            (top_scens / len(test_ref_df.dropna(subset=ref_wd_col)) * 100).to_frame(),
+            tablefmt="double_grid",
+            floatfmt=".1f",
+        ),
+    )
+
+    return test_ref_df.drop(columns=["rounded_wd", "all_turbines_waking"])