Wind Resource and Energy Yield Assessment: Hohn, Schleswig-Holstein

# --- Imports -------------------------------------------------------------

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from IPython.display import Markdown, display
from windforecast import config
from windforecast.io.dwd import load_hourly_wind

cfg = config.load()

pd.set_option("display.width", 120)
pd.set_option("display.max_columns", 25)

# --- Site map: turbine layout on satellite imagery ------------------------
import geopandas as gpd
import contextily as cx

from windforecast.paths import FIGURES

site = cfg["site"]

turbines = gpd.GeoDataFrame(
    {"name": [t["name"] for t in site["turbines"]]},
    geometry=gpd.points_from_xy(
        [t["longitude"] for t in site["turbines"]],
        [t["latitude"] for t in site["turbines"]],
    ),
    crs="EPSG:4326",
).to_crs(3857)

station = gpd.GeoDataFrame(
    {"name": [f"DWD {site['name']} (02303)"]},
    geometry=gpd.points_from_xy([site["longitude"]], [site["latitude"]]),
    crs="EPSG:4326",
).to_crs(3857)

fig, ax = plt.subplots(figsize=(10, 9))

# extent: bounding box of all points, padded
pad = 900  # m
xmin, ymin, xmax, ymax = gpd.pd.concat([turbines, station]).total_bounds
ax.set_xlim(xmin - pad, xmax + pad)
ax.set_ylim(ymin - pad, ymax + pad)

cx.add_basemap(ax, source=cx.providers.Esri.WorldImagery, crs=turbines.crs, zoom=14)

turbines.plot(ax=ax, marker="*", color="red", markersize=260,
              edgecolor="white", linewidth=1.0, label="Turbines T1-T5")
station.plot(ax=ax, marker="o", color="red", markersize=90, edgecolor="white",
             label="DWD station Hohn (02303)")

for pt, name in zip(turbines.geometry, turbines["name"]):
    ax.annotate(name, (pt.x, pt.y), xytext=(8, 6), textcoords="offset points",
                color="white", fontsize=10, fontweight="bold")

ax.annotate("Station", (station.geometry.iloc[0].x, station.geometry.iloc[0].y),
            xytext=(8, 6), textcoords="offset points",
            color="white", fontsize=10, fontweight="bold")

# title inside the frame
ax.set_title("")
ax.text(0.5, 0.975,
        f"Hypothetical five-turbine site at Hohn airfield\n"
        f"{site['turbine_count']} x {cfg['turbine']['model']}, "
        f"hub height {cfg['turbine']['hub_height_m']} m",
        transform=ax.transAxes, ha="center", va="top",
        color="white", fontsize=12, fontweight="bold", linespacing=1.4)

ax.set_xticks([])
ax.set_yticks([])
for spine in ax.spines.values():
    spine.set_visible(True)
    spine.set_color("white")
    spine.set_linewidth(1.2)

legend = ax.legend(loc="upper left", framealpha=0.85)

# scale bar, centered just below the legend box
bar_m = 1000
fig.canvas.draw()                                   # legend extent is known only after a draw
lb = legend.get_window_extent().transformed(ax.transAxes.inverted())
bar_frac = bar_m / (ax.get_xlim()[1] - ax.get_xlim()[0])
xc = (lb.x0 + lb.x1) / 2 - 0.03
xa0, xa1 = xc - bar_frac / 2, xc + bar_frac / 2
ya = lb.y0 - 0.025
ax.plot([xa0, xa1], [ya, ya], color="white", lw=3, transform=ax.transAxes)
ax.annotate("1 km", (xc, ya), xycoords="axes fraction", xytext=(0, 6),
            textcoords="offset points", ha="center", color="white", fontsize=9)

ax.annotate("Imagery: Esri World Imagery", (xmax + pad, ymin - pad),
            xytext=(-8, 8), textcoords="offset points", ha="right",
            color="white", fontsize=8, alpha=0.8)

fig.tight_layout()
fig.savefig(FIGURES / "site_map_hohn.png", dpi=500, bbox_inches="tight")
plt.show()

# --- Download DWD data:  -----------------------------------------------------

# from windforecast.download.dwd import main as download_dwd
# download_dwd()   # writes into data/raw/dwd/, skips files already there

# `load_hourly_wind` handles the DWD conventions: `-999` → NaN, MEZ → UTC, a complete hourly axis, 
# and direction codes `0` (calm) and `990` (variable) set to missing.

# Columns:
#   F      mean wind speed at 10 m (m/s)
#   D      wind direction, degrees (0/990 codes and calms -> NaN)
#   D_raw  direction as DWD delivers it
#   QN_3   DWD quality level of the measurement

dwd = cfg["reference_data"]["dwd"]
window = cfg["analysis"]["measurement_window"]
WINDOW = (str(window["start"]), str(window["end"]))

measured = load_hourly_wind(dwd["primary"],                       # load the station data
                            time_base=dwd["time_base"])           # convert MEZ to UTC
measured = measured.loc[WINDOW[0]:f"{WINDOW[1]} 23:00"]           # keep only the measurement window (2 years)


# check for missing data
hours_expected = len(pd.date_range(WINDOW[0], f"{WINDOW[1]} 23:00", freq="h"))
hours_speed = measured["F"].notna().sum()
hours_direction = measured["D"].notna().sum()

share_speed = 100 * hours_speed / hours_expected
share_direction = 100 * hours_direction / hours_expected
share_worst_month = 100 * (measured["F"].notna()
                           .groupby([measured.index.year, measured.index.month])
                           .mean().min())

display(Markdown(f"""
| | Hours w.<br>Measurements | Availability |
|---|---:|---:|
| Expected (2021–2022) | {hours_expected} | |
| With wind speed | {hours_speed} | {share_speed:.2f} % |
| With wind direction | {hours_direction} | {share_direction:.2f} % |
| Worst calendar month | | {share_worst_month:.2f} % |

Availability is {share_speed:.2f} % for speed and {share_direction:.2f} % for direction, comfortably above 80 %. So our "campaign" went very well!


"""))

from windforecast.analysis.profile import add_z0, log_extrapolate

z0_sectors = [(s["from_deg"], s["to_deg"], s["z0_m"]) for s in cfg["analysis"]["roughness"]["sectors"]]
measured = log_extrapolate(df=add_z0(df=measured, sectors=z0_sectors),
                           speed_col="F", z0_col="z0", from_height=10, to_height=100)

display(measured[["F", "F100_extrapol"]].describe())

u10, u100 = measured["F"].mean(), measured["F100_extrapol"].mean()
display(Markdown(f'We have our measured wind at 100 m now. The mean wind at 10 m is **{u10:.2f} m/s**, and at 100 m it is **{u100:.2f} m/s**, showing why the higher the better for wind turbines. '
                 f'Two years is nothing in climate terms, though, and one windy year can make a site look great! But the turbines will run for about 20 years, meaning we need a long-term record to put our short measurement into context. '
                 f'And that is exactly what Steps 2 and 3 do!'))

# --- Download COSMO SP_100 and DD_100m: -----------------------------------------
# from windforecast.download.cosmo import download_group

# download_group(
#     cfg["reference_data"]["cosmo"]["base_url"],
#     "hourly_zlev",
#     ["SP_100m", "DD_100m"],
#     site["latitude"], site["longitude"],
# )   # skips months already in the CSV; writes data/raw/cosmo/hourly_zlev/<var>/cosmo_<var>.csv

from windforecast.paths import RAW_COSMO_HOURLY

sp = pd.read_csv(RAW_COSMO_HOURLY / "SP_100m" / "cosmo_SP_100m.csv", index_col="time", parse_dates=True)
dd = pd.read_csv(RAW_COSMO_HOURLY / "DD_100m" / "cosmo_DD_100m.csv", index_col="time", parse_dates=True)

cosmo = pd.concat([sp.iloc[:, 0].rename("SP_100m"),
                   dd.iloc[:, 0].rename("DD_100m")], axis=1)

# 12 sectors of 30°, sector 0 = 345°–15°
cosmo["sector"] = (((cosmo["DD_100m"] + 15) % 360) // 30).astype("Int64")
display(cosmo.describe())

u_cosmo_lt = cosmo["SP_100m"].mean()
u_cosmo_window = cosmo.loc[WINDOW[0]:f"{WINDOW[1]} 23:00", "SP_100m"].mean()

display(Markdown(f"""
Now our COSMO data is ready. It gives us {len(cosmo):,} hours of wind at 100 m, from {cosmo.index.min():%Y} to {cosmo.index.max():%Y},
with a mean of **{u_cosmo_lt:.2f} m/s**. Over our two measurement years alone it averages {u_cosmo_window:.2f} m/s,
about {100 * (u_cosmo_window / u100 - 1):.0f} % more than our extrapolated measurement ({u100:.2f} m/s) over the same period.
For the MCP itself the gap doesn't matter: the regression is fitted to our measurement, so COSMO only tells us how windy the long term was *compared to* 2021–2022.
But a gap this size is relevant and also a hint that the extrapolation may underestimate the wind. More on that in Step 4.
"""))

# Joining measurement and COSMO in one data frame:

short = measured["F100_extrapol"]           # measurement
ref = cosmo["SP_100m"]                      # COSMO
sect = cosmo["sector"]

both = pd.concat({"T": short, "R": ref, "sect": sect}, axis=1).dropna()

import statsmodels.api as sm

# correlate: one linear regression per sector
rows = {}
for s, g in both.groupby("sect"):
    model = sm.OLS(g["T"], sm.add_constant(g["R"])).fit()
    rows[s] = {"a": model.params["const"], "b": model.params["R"],
               "r": g["T"].corr(g["R"]), "n": len(g)}
coef = pd.DataFrame(rows).T

display(coef.round(3))

# predict: each long-term hour gets the a, b of its sector
cosmo_sec = cosmo["sector"]
predicted_sec = pd.DataFrame({
    "u": (cosmo_sec.map(coef["a"]) + cosmo_sec.map(coef["b"]) * cosmo["SP_100m"]).clip(lower=0),
    "sector": cosmo_sec,
})


# numbers for the text
NAMES = ["north", "north-north-east", "east-north-east", "east", "east-south-east", "south-south-east",
         "south", "south-south-west", "west-south-west", "west", "west-north-west", "north-north-west"]
s_min, s_max = int(coef["r"].idxmin()), int(coef["r"].idxmax())
r_w = (coef["r"] * coef["n"]).sum() / coef["n"].sum()
share_good = 100 * coef.loc[coef["r"] >= 0.8, "n"].sum() / coef["n"].sum()
u_lt = predicted_sec["u"].mean()
diff = 100 * (u_lt / u100 - 1)

display(Markdown(f"""
Per sector, r ranges from {coef['r'].min():.2f} (sector {s_min}, {NAMES[s_min]}) to {coef['r'].max():.2f}
(sector {s_max}, {NAMES[s_max]}). Weighted by the number of hours, the mean r is {r_w:.2f},
{'a bit below' if r_w < 0.8 else 'above'} the 0.8 that MEASNET 8.3.2 recommends, and {share_good:.0f} % of the hours
fall in sectors that reach it. Not surprising for a 10 m station lifted to 100 m with a formula,
but worth remembering for the uncertainty budget later.

So, were our two years normal? The measurement averaged {u100:.2f} m/s at 100 m, and the 24-year long-term series
averages {u_lt:.2f} m/s. The measurement period was about {abs(diff):.0f} % {'calmer' if diff > 0 else 'windier'}
than usual, so without this step we would have {'under' if diff > 0 else 'over'}estimated the wind by about that much,
and the energy by roughly twice that.

Here's what the long-term wind looks like: where it comes from and how strong it is.
"""))

from windforecast.analysis.profile import windrose, speed_distribution

# --- both panels, one figure -------------------------------------------------
SPEED_BINS = (0, 4, 6, 8, 10, 12, np.inf)

fig = plt.figure(figsize=(14, 6))
ax_rose = fig.add_subplot(1, 2, 1, projection="polar")
ax_hist = fig.add_subplot(1, 2, 2)

_, rose = windrose(cosmo["sector"], predicted_sec["u"], speed_bins=SPEED_BINS,
                   ax=ax_rose, title="Richtungsverteilung")
_, dist = speed_distribution(predicted_sec["u"], bin_width=0.25, speed_bins=SPEED_BINS,
                             ax=ax_hist, title="Geschwindigkeitsverteilung")

# the rose already carries the class legend; keep one for the figure
ax_rose.legend(title="u [m/s]", loc="upper left", bbox_to_anchor=(-0.28, 1.05),
               frameon=False, fontsize=9)

fig.suptitle("Hohn, 100 m — langzeitkorrigiert (COSMO-R6G2)", fontsize=13)
fig.tight_layout()
plt.show()

freq = 100 * predicted_sec["sector"].value_counts(normalize=True).sort_index()
display(Markdown(f"""
{freq.loc[7:10].sum():.0f} % of all hours come from between south-south-west and west-north-west (195°–315°), and that's
also where the strongest winds are, with sector 10 (WNW) averaging {predicted_sec.groupby('sector')['u'].mean()[10]:.1f} m/s.
There's a smaller second peak from the east-north-east.
"""))

# --- Download GWA wind speed: -----------------------------------------------
# from windforecast.download.gwa import main as download_gwa

# download_gwa()   # writes DEU_wind-speed_<h>m.tif into data/raw/gwa/, skips what is there

from windforecast.io.gwa import transfer_factors

TARGET_HEIGHT = 100   # hub height

turbines = [(t["latitude"], t["longitude"]) for t in site["turbines"]]
names = [t["name"] for t in site["turbines"]]

transfer = transfer_factors(
    (site["latitude"], site["longitude"]), turbines,
    height=TARGET_HEIGHT, names=names,
)
display(transfer)

# numbers for the text

max = transfer["factor"].idxmax()
min = transfer["factor"].idxmin()
min_val = 100 * (transfer.loc[min, "factor"] - 1)
max_val = 100 * (transfer.loc[max, "factor"] - 1)
gwa_station = transfer["gwa_station_100m"].iloc[0]
u_lt = predicted_sec["u"].mean()
gap = 100 * (1 - u_lt / gwa_station)

display(Markdown(f"""
Each turbine now has a factor: how much stronger or weaker the wind is there than at the station, according to GWA.
They're all within about ±1 % ({min} {min_val:+.1f} %, {max} {max_val:+.1f} %), which makes sense over open, flat grassland.
Small, but with the energy factor of 1.5–2.5 from above, ±1 % in wind is about ±2 % in energy.
These factors come back in Step 6.

A note on levels: GWA says {gwa_station:.2f} m/s at 100 m at the station, and our long-term series says
{u_lt:.2f} m/s, about {gap:.0f} % lower. Part of that is simply different models and periods (GWA covers 2008–2017),
and it's exactly why I only take ratios from GWA. But it's also a warning about Step 1. Lázár et al. [5] tested a very
similar log-profile extrapolation against measurements and found it underestimated the wind more than 70 % of the time,
with the largest misses at night under stable conditions, when, as I read it, the wind higher up decouples from the ground.
They extrapolated from 50 m; from 10 m the effect is probably larger, since the extrapolation uncertainty grows with the
height ratio (DFBEW [1] Table 2: 1.1–6.3 %; Waars [4] §2.1.4). Their errors are smaller than our gap, so this points in
the same direction rather than explaining all of it. Either way, read the absolute energy numbers at the end as conservative.

DFBEW [1] (III.1) is clear about this kind of error: *"Systematische Fehler dürfen bei der Windpotenzialbestimmung
nicht auftreten oder müssen korrigiert werden."* Systematic errors must not occur, or they have to be corrected.
MEASNET [2] 8.6 agrees: a known bias belongs in its own correction, not in the uncertainty. I haven't corrected it yet,
but I'll come back to it in a later update of this notebook.
"""))

# --- Download COSMO 1hrPt_tas, 1hrPt_ps and 1hrPt_huss: -------------------------------------
# download_group(
#     cfg["reference_data"]["cosmo"]["base_url"],
#     "thermo",
#     ["1hrPt_tas", "1hrPt_ps", "1hrPt_huss"],
#     site["latitude"], site["longitude"],
# )   # skips months already in the CSV; writes data/raw/cosmo/hourly/<var>/cosmo_<var>.csv

from windforecast.paths import RAW_COSMO_THERMO

parts = []
for var in ["1hrPt_tas", "1hrPt_ps", "1hrPt_huss"]:
    df = pd.read_csv(RAW_COSMO_THERMO / var / f"cosmo_{var}.csv", index_col="time", parse_dates=True)
    parts.append(df.iloc[:, 0].rename(var.replace("1hrPt_", "")))

cosmo_thermo = pd.concat(parts, axis=1).sort_index()

print(f"{cosmo_thermo.index.min()} .. {cosmo_thermo.index.max()}  "
      f"({len(cosmo_thermo)} rows)")

expected = pd.date_range(cosmo_thermo.index.min(), cosmo_thermo.index.max(), freq="h")
print(f"expected hours {len(expected)}, missing timestamps {len(expected.difference(cosmo_thermo.index))}")
print(f"duplicated timestamps {cosmo_thermo.index.duplicated().sum()}\n")

print("availability per column (% of rows with a value):")
print((100 * cosmo_thermo.notna().mean()).round(2).to_string())
print()
print(cosmo_thermo.describe().round(4).to_string())


display(Markdown(f"""
Everything here looks good!
"""))

R_D, G0, L = 287.05, 9.80665, -0.0065

def density_at(df, from_height, to_height):
    Tv1 = df["tas"] * (1 + 0.61 * df["huss"])
    dz = to_height - from_height
    Tv2 = Tv1 + L * dz
    P2 = df["ps"] * (Tv2 / Tv1) ** (-G0 / (R_D * L))
    return P2 / (R_D * Tv2)

cosmo_thermo["rho_100m"] = density_at(cosmo_thermo, from_height=2.0, to_height=100.0)
predicted_sec["rho_100m"] = cosmo_thermo["rho_100m"]
rho_mean = predicted_sec["rho_100m"].mean()

display(predicted_sec.head())

display(Markdown(f"""
So this is what our table looks like. On average, the air at 100 m is **{rho_mean:.3f} kg/m³**.
The next step is to plug it into the equation for each turbine (and don't forget the turbine factor from Step 4)
to get the energy yield of our wind farm!
"""))

# fit one Weibull distribution per sector to the predicted long-term wind
from scipy import stats

rows = {}
for s, g in predicted_sec.groupby("sector"):
    u = g["u"][g["u"] > 0]
    k, _, A = stats.weibull_min.fit(u, floc=0)
    rows[s] = {"A": A, "k": k, "n": len(g),
               "freq": len(g) / len(predicted_sec),
               "u_mean": g["u"].mean(),
               "rho": g["rho_100m"].mean()}

weibull = pd.DataFrame(rows).T
display(weibull.round(4))


display(Markdown(f"""
This table is the input for the energy calculation. For each of the {len(weibull)} sectors it holds the Weibull
parameters (A, k), how often the wind comes from that sector (freq) and the mean air density there (rho),
all based on {int(weibull['n'].sum()):,} long-term hours."""))

# Nordex N100/2500 power curve at 1.225 kg/m³ (Staffell (no year) [6])
from windforecast.paths import CONFIG

pc = pd.read_csv(CONFIG / cfg["turbine"]["power_curve"])
ws, kw = pc["wind_speed_ms"].to_numpy(), pc["power_kw"].to_numpy()

RHO_REF = cfg["turbine"]["power_curve_reference_density"]   # density the curve is valid for
P_RATED = kw.max()                                          # rated power in kW
CUT_OUT = ws[kw > 0].max()                                  # 25 m/s

def power_curve(u):
    u = np.asarray(u, dtype=float)
    return np.where(u > CUT_OUT, 0.0, np.interp(u, ws, kw, left=0.0, right=0.0))

display(Markdown(f"""
Here we read the power curve and create a function that takes any set of wind speeds and returns the power the
turbine would produce at each, interpolating linearly between the tabulated points. 
"""))

bins = np.arange(0, 30.5, 1.0)                       # 1 m/s wind speed bins
centres = (bins[:-1] + bins[1:]) / 2

def turbine_aep(factor=1.0):
    """Gross AEP in kWh/yr from the sector Weibull fits."""
    aep = 0.0
    for s, row in weibull.iterrows():
        cdf = stats.weibull_min.cdf(bins, row["k"], loc=0, scale=row["A"] * factor)
        p = np.diff(cdf)                                         # probability of each bin
        u_eq = centres * (row["rho"] / RHO_REF) ** (1/3)         # density correction (Step 5)
        aep += row["freq"] * (p * power_curve(u_eq)).sum()       # weighted by sector frequency
    return aep * 8766

# all five turbines, each with its own GWA factor (Step 4)
transfer["aep_MWh"] = [turbine_aep(f) / 1e3 for f in transfer["factor"]]
farm_GWh = transfer["aep_MWh"].sum() / 1e3
spread = 100 * (transfer.loc["T1", "aep_MWh"] / transfer.loc["T5", "aep_MWh"] - 1)

# one turbine at the station, for comparison
aep_station = turbine_aep(1.0)
flh = aep_station / P_RATED

display(Markdown(f"""
MEASNET 10.2: the gross yield is calculated per turbine, and the farm is simply the sum. This is where the GWA factors
from Step 4 come in. Multiplying every wind speed by a factor f stretches the Weibull distribution: A becomes f·A,
and k stays the same. So each turbine gets the same calculation with its own A.
"""))

display(transfer[["distance_km", "factor", "aep_MWh"]].round({"distance_km": 2, "factor": 3, "aep_MWh": 1}))

display(Markdown(f"""
The farm produces **{farm_GWh:.1f} GWh per year** gross. T1 is the best spot and T5 the weakest,
a spread of {spread:.1f} %.

For comparison, one turbine standing right at the station would produce {aep_station / 1e6:.2f} GWh per year.
That's about {flh:,.0f} full-load hours: the number of hours it would have to run at its full rated power
({P_RATED:,.0f} kW) to produce the same energy. Divided by the 8766 hours in a year, that gives the capacity factor
of {100 * flh / 8766:.0f} %, the share of its maximum possible output the turbine actually delivers.
"""))

# time-series check (MEASNET 10.2): power curve straight over the long-term hourly series
u_ts = predicted_sec["u"]
rho_ts = predicted_sec["rho_100m"].fillna(predicted_sec["rho_100m"].mean())   # 2 missing hours
years = len(predicted_sec) / 8766

def turbine_aep_ts(factor=1.0):
    """Gross AEP in kWh/yr from the hourly long-term series."""
    return power_curve(factor * u_ts * (rho_ts / RHO_REF) ** (1/3)).sum() / years

farm_ts_GWh = sum(turbine_aep_ts(f) for f in transfer["factor"]) / 1e6
check = 100 * (farm_GWh / farm_ts_GWh - 1)

display(Markdown(f"""
The time-series check gives {farm_ts_GWh:.1f} GWh per year for the farm, a difference of {check:+.1f} %,
{'small enough to trust the Weibull fit for the energy calculation' if abs(check) < 2 else 'large enough to look into the Weibull fit before going further'}.
"""))