"""Curated, source-backed European offshore reference farms."""
from __future__ import annotations
import copy
import hashlib
import json
import os
import tempfile
import time
import unicodedata
import urllib.request
from dataclasses import dataclass, replace
from functools import lru_cache
from pathlib import Path
from typing import Literal
import numpy as np
import shapely
from ._cluster import Cluster
from ._validation import _validate_integer
from .boundaries.base import _GeometryBoundary
from .layouts.base import Farm
from .resources import GWA
from .turbines import Tabulated, Turbine
__all__ = ["ReferenceFarmInfo", "ReferenceCatalog", "reference_farm"]
_CATALOG_PATH = Path(__file__).parent / "data" / "reference_farms.npz"
_CITATIONS = (
"Fischereit, Vollmer & Hansen, Open European offshore wind turbine "
"database, doi:10.5281/zenodo.19879819",
"OpenStreetMap contributors, current and historical turbine nodes, ODbL 1.0",
"EMODnet Human Activities wind-farm polygons, CC BY 4.0",
)
_GWA_DOWNLOAD_TIMEOUT = 30.0
_GWA_MAX_PAYLOAD_BYTES = 1024 * 1024
[docs]
@dataclass(frozen=True)
class ReferenceFarmInfo:
"""Searchable metadata for one packaged farm or named farm group.
Coordinates of the centroid are WGS 84 degrees. ``boundary_kind`` is
``"area"``, ``"linear"`` for a single row of turbines, or ``"cluster"`` for
a group. ``boundary_containment`` is the share of turbines inside the
packaged boundary, ``boundary_containment_50m`` the share within 50 m of
it and ``boundary_max_outside_distance_m`` the farthest turbine outside it,
in meters. ``quality_flags`` name known caveats such as positions outside
the boundary. ``fit_eligible`` marks an area farm of at least five distinct
positions with 95 percent of them within 50 m of the boundary, and
``shape_fit_eligible`` an area farm whose boundary is one valid EMODnet
polygon. A group reports its members' weakest containment and is eligible
for neither fit.
"""
name: str
aliases: tuple[str, ...]
country: str
commissioning_year: int | None
latitude: float
longitude: float
n_turbines: int
installation_type: str
boundary_source: str
boundary_kind: str
boundary_containment: float
boundary_containment_50m: float
boundary_max_outside_distance_m: float
boundary_source_names: tuple[str, ...]
boundary_source_feature_ids: tuple[str, ...]
boundary_adjustment: dict[str, object]
boundary_vertex_count_raw: int
quality_flags: tuple[str, ...]
fit_eligible: bool
shape_fit_eligible: bool
source_dates: dict[str, str]
is_group: bool = False
members: tuple[str, ...] = ()
catalog_date: str = ""
citations: tuple[str, ...] = _CITATIONS
@lru_cache(maxsize=1)
def _catalog() -> dict[str, np.ndarray]:
if not _CATALOG_PATH.exists():
raise RuntimeError(
"OOWL's reference-farm catalog is missing; reinstall OOWL with "
"its package data"
)
with np.load(_CATALOG_PATH, allow_pickle=False) as data:
arrays = {key: data[key] for key in data.files}
if int(arrays.get("schema_version", 0)) != 5:
raise RuntimeError("unsupported reference-farm catalog schema")
return arrays
def _json_tuple(value) -> tuple[str, ...]:
return tuple(str(item) for item in json.loads(str(value)))
def _year(value) -> int | None:
"""Commissioning year, stored as zero or less when unknown."""
year = int(value)
return None if year <= 0 else year
def _member_indices(data, members: tuple[str, ...]) -> list[int]:
farm_names = list(map(str, data["farm_name"]))
return [farm_names.index(name) for name in members]
def _atomic_info(index: int) -> ReferenceFarmInfo:
data = _catalog()
return ReferenceFarmInfo(
name=str(data["farm_name"][index]),
aliases=_json_tuple(data["aliases_json"][index]),
country=str(data["country"][index]),
commissioning_year=_year(data["commissioning_year"][index]),
latitude=float(data["centroid_latlon"][index, 0]),
longitude=float(data["centroid_latlon"][index, 1]),
n_turbines=int(data["n_turbines"][index]),
installation_type=str(data["installation_type"][index]),
boundary_source=str(data["boundary_source"][index]),
boundary_kind=str(data["boundary_kind"][index]),
boundary_containment=float(data["boundary_containment"][index]),
boundary_containment_50m=float(data["boundary_containment_50m"][index]),
boundary_max_outside_distance_m=float(
data["boundary_max_outside_distance_m"][index]
),
boundary_source_names=_json_tuple(data["boundary_source_names_json"][index]),
boundary_source_feature_ids=_json_tuple(
data["boundary_source_feature_ids_json"][index]
),
boundary_adjustment=json.loads(str(data["boundary_adjustment_json"][index])),
boundary_vertex_count_raw=int(data["boundary_vertex_count_raw"][index]),
quality_flags=_json_tuple(data["quality_flags_json"][index]),
fit_eligible=bool(data["fit_eligible"][index]),
shape_fit_eligible=bool(data["shape_fit_eligible"][index]),
source_dates={
"catalog" if key == "catalogue" else str(key): str(value)
for key, value in json.loads(str(data["source_dates_json"][index])).items()
},
catalog_date=str(data["catalogue_date"]),
)
def _group_info(index: int) -> ReferenceFarmInfo:
"""Group metadata, aggregated from its member farms."""
data = _catalog()
members = _json_tuple(data["group_members_json"][index])
member_indices = _member_indices(data, members)
def union(key: str) -> tuple[str, ...]:
"""Members' values for ``key`` in first-seen order, without repeats."""
return tuple(
dict.fromkeys(
value
for member in member_indices
for value in _json_tuple(data[key][member])
)
)
member_source_dates = [
json.loads(str(data["source_dates_json"][member])) for member in member_indices
]
emodnet_years = sorted(
{
year
for source_dates in member_source_dates
for year in source_dates["emodnet_update_years"].split(", ")
if year
}
)
osm_dates = [
source_dates["osm_latest_selected_object"]
for source_dates in member_source_dates
if source_dates["osm_latest_selected_object"]
]
return ReferenceFarmInfo(
name=str(data["group_name"][index]),
aliases=_json_tuple(data["group_aliases_json"][index]),
country=str(data["group_country"][index]),
commissioning_year=_year(data["group_commissioning_year"][index]),
latitude=float(data["group_centroid_latlon"][index, 0]),
longitude=float(data["group_centroid_latlon"][index, 1]),
n_turbines=int(data["group_n_turbines"][index]),
installation_type="Grounded",
boundary_source="EMODnet",
boundary_kind="cluster",
boundary_containment=float(
np.min(data["boundary_containment"][member_indices])
),
boundary_containment_50m=float(
np.min(data["boundary_containment_50m"][member_indices])
),
boundary_max_outside_distance_m=float(
np.max(data["boundary_max_outside_distance_m"][member_indices])
),
boundary_source_names=union("boundary_source_names_json"),
boundary_source_feature_ids=union("boundary_source_feature_ids_json"),
boundary_adjustment={"kind": "none", "source_geometry_preserved": True},
boundary_vertex_count_raw=int(
np.sum(data["boundary_vertex_count_raw"][member_indices])
),
quality_flags=("exact_member_boundaries",),
fit_eligible=False,
shape_fit_eligible=False,
source_dates={
"catalog": str(data["catalogue_date"]),
"eww_cohort_cutoff": "2026-01-27",
"emodnet_update_years": ", ".join(emodnet_years),
"osm_latest_selected_object": max(osm_dates, default=""),
},
is_group=True,
members=members,
catalog_date=str(data["catalogue_date"]),
)
def _normalize_name(value: str) -> str:
"""Casefolded letters and digits only, with accents removed."""
text = unicodedata.normalize("NFKD", str(value))
text = "".join(char for char in text if not unicodedata.combining(char))
return "".join(char.casefold() for char in text if char.isalnum())
@lru_cache(maxsize=1)
def _name_index() -> tuple[tuple[str, int, str, frozenset[str]], ...]:
"""``(kind, index, canonical, normalized names)`` for every farm and group.
Cached once since the packaged catalog never changes.
"""
data = _catalog()
return tuple(
(
kind,
index,
str(canonical),
frozenset(
_normalize_name(candidate)
for candidate in (str(canonical), *_json_tuple(alias_json))
),
)
for kind, names, aliases in (
("farm", data["farm_name"], data["aliases_json"]),
("group", data["group_name"], data["group_aliases_json"]),
)
for index, (canonical, alias_json) in enumerate(zip(names, aliases))
)
@lru_cache(maxsize=2)
def _catalog_infos(include_groups: bool) -> tuple[ReferenceFarmInfo, ...]:
"""Metadata of every farm, then every group if ``include_groups``.
Cached since the packaged catalog never changes. Callers copy the
dictionary fields before handing an entry out.
"""
data = _catalog()
infos = tuple(_atomic_info(index) for index in range(len(data["farm_name"])))
if include_groups:
infos += tuple(_group_info(index) for index in range(len(data["group_name"])))
return infos
def _resolve_name(name: str) -> tuple[str, int]:
"""Return ``("farm" | "group", index)`` for an exact name or alias."""
query = _normalize_name(name)
if not query:
raise ValueError("farm name must not be empty")
exact = []
partial = []
for kind, index, canonical, normalized in _name_index():
if query in normalized:
exact.append((kind, index, canonical))
elif any(query in item or item in query for item in normalized):
partial.append(canonical)
if len(exact) == 1:
kind, index, _canonical = exact[0]
return kind, index
if len(exact) > 1:
matches = sorted(canonical for _kind, _index, canonical in exact)
raise ValueError(
f"ambiguous reference farm {name!r}; matches: {', '.join(matches)}"
)
matches = sorted(set(partial))
if not matches:
raise KeyError(f"unknown reference farm {name!r}")
raise ValueError(
f"reference farm {name!r} is not an exact name or alias; "
f"possible matches: {', '.join(matches)}"
)
def _is_local(coordinate_frame: str) -> bool:
if coordinate_frame == "local":
return True
if coordinate_frame.lower() in {"epsg:3035", "3035"}:
return False
raise ValueError("coordinate_frame must be 'local' or 'epsg:3035'")
def _slice(offsets, index):
return slice(int(offsets[index]), int(offsets[index + 1]))
def _wkb_geometry(data, prefix: str, index: int):
"""Boundary ``index`` of the packed WKB arrays ``prefix``, or ``None``."""
segment = _slice(data[f"{prefix}_wkb_offsets"], index)
payload = data[f"{prefix}_wkb_data"][segment].tobytes()
if not payload:
return None
geometry = shapely.from_wkb(payload)
if geometry is None or geometry.is_empty:
raise ValueError("reference boundary is empty")
return geometry
def _turbine_types(data, type_indices: np.ndarray):
"""One shared turbine, or one per position when the types differ."""
types = {}
for key in sorted(set(type_indices.tolist())):
curve = _slice(data["curve_offsets"], key)
types[key] = Tabulated(
rotor_diameter=float(data["type_rotor_diameter"][key]),
hub_height=float(data["type_hub_height"][key]),
wind_speeds=data["curve_wind_speed"][curve],
power_values=data["curve_power_w"][curve],
ct_values=data["curve_ct"][curve],
name=str(data["type_name"][key]),
rated_power=float(data["type_rated_power_w"][key]),
)
assigned = tuple(types[key] for key in type_indices.tolist())
if assigned and all(item is assigned[0] for item in assigned):
return assigned[0]
return assigned
def _optional_int(value) -> int | None:
"""OpenStreetMap identifiers, stored as negative when absent."""
value = int(value)
return None if value < 0 else value
def _position_records(data, indices) -> list[dict]:
return [
{
"canonical_farm": str(data["position_canonical_farm"][index]),
"longitude": float(data["positions_wgs84"][index, 0]),
"latitude": float(data["positions_wgs84"][index, 1]),
"position_source": str(data["position_source"][index]),
"position_source_id": str(data["position_source_id"][index]),
"osm_id": _optional_int(data["position_osm_id"][index]),
"osm_version": _optional_int(data["position_osm_version"][index]),
"osm_timestamp": str(data["position_osm_timestamp"][index]),
"osm_valid_to": str(data["position_osm_valid_to"][index]),
"reference": str(data["position_reference"][index]),
"lifecycle_status": str(data["position_lifecycle_status"][index]),
"status_source": str(data["position_status_source"][index]),
"status_date": str(data["position_status_date"][index]),
"status_evidence": str(data["position_status_evidence"][index]),
}
for index in map(int, indices)
]
def _cache_path(
info: ReferenceFarmInfo, roughness_length: float, reference_height: float
) -> Path:
key = (
f"{info.latitude:.6f}|{info.longitude:.6f}|"
f"{float(roughness_length):.8g}|{float(reference_height):.4f}|gwa-v4"
)
digest = hashlib.sha256(key.encode()).hexdigest()[:20]
# The XDG spec treats an empty value as unset.
root = Path(os.environ.get("XDG_CACHE_HOME") or Path.home() / ".cache")
return root / "oowl" / "gwa" / f"{digest}.lib"
def _cached_gwa(
path: Path,
info: ReferenceFarmInfo,
roughness_length: float,
reference_height: float,
) -> GWA:
return GWA(
info.latitude,
info.longitude,
roughness_length,
reference_height,
libfile_path=path,
)
def _download_gwa(
cache: Path,
info: ReferenceFarmInfo,
roughness_length: float,
reference_height: float,
) -> GWA:
"""Download the libfile, validate it, then publish it as ``cache``.
The download goes to a private temporary file. Linking it into place never
replaces an entry that a concurrent download published first.
"""
url = f"{GWA.GWA_URL}?lat={info.latitude}&long={info.longitude}&areaname=oowl"
request = urllib.request.Request(url, headers=GWA.REQUEST_HEADERS)
cache.parent.mkdir(parents=True, exist_ok=True)
temporary = tempfile.NamedTemporaryFile(
mode="wb",
dir=cache.parent,
prefix=f".{cache.name}.",
suffix=".tmp",
delete=False,
)
temporary_path = Path(temporary.name)
try:
deadline = time.monotonic() + _GWA_DOWNLOAD_TIMEOUT
with (
temporary,
urllib.request.urlopen(
request, context=GWA._ssl_context(), timeout=_GWA_DOWNLOAD_TIMEOUT
) as response,
):
content_length = response.headers.get("Content-Length")
if (
content_length is not None
and int(content_length) > _GWA_MAX_PAYLOAD_BYTES
):
raise RuntimeError("Global Wind Atlas response is too large")
downloaded = 0
while chunk := response.read(64 * 1024):
downloaded += len(chunk)
if downloaded > _GWA_MAX_PAYLOAD_BYTES:
raise RuntimeError("Global Wind Atlas response is too large")
if time.monotonic() > deadline:
raise TimeoutError("Global Wind Atlas download timed out")
temporary.write(chunk)
resource = _cached_gwa(temporary_path, info, roughness_length, reference_height)
try:
os.link(temporary_path, cache)
except FileExistsError:
return _cached_gwa(cache, info, roughness_length, reference_height)
resource.libfile_path = cache
return resource
finally:
temporary_path.unlink(missing_ok=True)
def _gwa_resource(
info: ReferenceFarmInfo, roughness_length: float, reference_height: float
) -> GWA:
"""Global Wind Atlas resource at the farm centroid, cached on disk.
An unreadable cache entry is removed and downloaded again.
"""
cache = _cache_path(info, roughness_length, reference_height)
try:
return _cached_gwa(cache, info, roughness_length, reference_height)
except Exception:
cache.unlink(missing_ok=True)
return _download_gwa(cache, info, roughness_length, reference_height)
def _reference_metadata(info, coordinate_frame, crs, indices):
return {
"kind": "reference_farm",
"name": info.name,
"layout_basis": "all_commissioned_positions",
"country": info.country,
"commissioning_year": info.commissioning_year,
"catalog_date": info.catalog_date,
"centroid_latitude": info.latitude,
"centroid_longitude": info.longitude,
"coordinate_frame": coordinate_frame,
"crs": crs,
"installation_type": info.installation_type,
"boundary_source": info.boundary_source,
"boundary_source_names": list(info.boundary_source_names),
"boundary_source_feature_ids": list(info.boundary_source_feature_ids),
"boundary_adjustment": dict(info.boundary_adjustment),
"boundary_kind": info.boundary_kind,
"boundary_containment": info.boundary_containment,
"boundary_containment_50m": info.boundary_containment_50m,
"boundary_max_outside_distance_m": info.boundary_max_outside_distance_m,
"quality_flags": list(info.quality_flags),
"source_dates": dict(info.source_dates),
"turbine_records": _position_records(_catalog(), indices),
"citations": list(info.citations),
"modifiers": [],
}
def _reference_farm(
info,
positions,
geometry,
crs,
indices,
coordinate_frame,
resource,
roughness_length,
reference_height,
):
"""Build one farm and attach the requested resource."""
data = _catalog()
turbines = _turbine_types(data, data["position_type_index"][indices].astype(int))
boundary = None if geometry is None else _GeometryBoundary(geometry)
farm = Farm(
positions,
boundary,
metadata=_reference_metadata(info, coordinate_frame, crs, indices),
turbines=turbines,
)
if farm.n_turbines != info.n_turbines:
raise RuntimeError("reference catalog count does not match loaded positions")
if resource is None:
return farm
if reference_height is None:
assigned = (turbines,) if isinstance(turbines, Turbine) else turbines
reference_height = np.median([item.hub_height for item in assigned])
farm.resource = _gwa_resource(
info, float(roughness_length), float(reference_height)
)
return farm
def _load_atomic(
index, coordinate_frame, resource, roughness_length, reference_height
) -> Farm:
data = _catalog()
indices = np.arange(
int(data["position_offsets"][index]),
int(data["position_offsets"][index + 1]),
dtype=np.int64,
)
if _is_local(coordinate_frame):
positions = data["positions_local"][indices].astype(float, copy=True)
geometry = _wkb_geometry(data, "boundary_local", index)
crs = str(data["local_crs"][index])
else:
positions = data["positions_epsg3035"][indices].astype(float, copy=True)
geometry = _wkb_geometry(data, "boundary_epsg3035", index)
crs = "EPSG:3035"
return _reference_farm(
_atomic_info(index),
positions,
geometry,
crs,
indices,
coordinate_frame,
resource,
roughness_length,
reference_height,
)
def _load_group(
index, coordinate_frame, resource, roughness_length, reference_height
) -> Cluster:
data = _catalog()
group = _group_info(index)
group_slice = _slice(data["group_position_offsets"], index)
position_indices = data["group_position_indices"][group_slice].astype(np.int64)
# Local member boundaries are packed group after group.
first_member = sum(
len(_json_tuple(value)) for value in data["group_members_json"][:index]
)
# Members share the group's local frame, so their positions stay comparable.
local = _is_local(coordinate_frame)
if local:
all_positions = data["group_positions_local"][group_slice].astype(
float, copy=True
)
crs = str(data["group_local_crs"][index])
else:
all_positions = data["positions_epsg3035"][position_indices].astype(
float, copy=True
)
crs = "EPSG:3035"
farms = []
start = 0
for member, farm_index in enumerate(_member_indices(data, group.members)):
info = _atomic_info(farm_index)
stop = start + info.n_turbines
if local:
geometry = _wkb_geometry(
data, "group_member_boundary_local", first_member + member
)
else:
geometry = _wkb_geometry(data, "boundary_epsg3035", farm_index)
farms.append(
_reference_farm(
info,
all_positions[start:stop],
geometry,
crs,
position_indices[start:stop],
coordinate_frame,
resource,
roughness_length,
reference_height,
)
)
start = stop
cluster = Cluster(
tuple(farms),
metadata={
"kind": "reference_cluster",
"name": group.name,
"members": list(group.members),
"country": group.country,
"commissioning_year": group.commissioning_year,
"catalog_date": group.catalog_date,
"centroid_latitude": group.latitude,
"centroid_longitude": group.longitude,
"coordinate_frame": coordinate_frame,
"crs": crs,
"citations": list(group.citations),
},
)
if cluster.n_turbines != group.n_turbines:
raise RuntimeError("reference catalog count does not match loaded positions")
return cluster
[docs]
class ReferenceCatalog:
"""Callable access to packaged offshore farms and catalog metadata."""
[docs]
def __call__(
self,
name: str,
*,
coordinate_frame: str = "local",
resource: Literal["gwa"] | None = None,
roughness_length: float = 0.0002,
reference_height: float | None = None,
) -> Farm | Cluster:
"""Load one real farm, or a group as a :class:`~oowl.Cluster` of its members.
Every commissioned, source-backed turbine position is included.
Parameters
----------
name
Farm or group name or alias, ignoring case, accents, spaces and
punctuation.
coordinate_frame
``"local"`` for the farm's own projection in meters, or
``"epsg:3035"`` for the ETRS89 Lambert equal-area grid of Europe,
also in meters.
resource
``None`` loads offline. ``"gwa"`` downloads the Global Wind Atlas
climate at each farm centroid and caches it under
``$XDG_CACHE_HOME/oowl/gwa``.
roughness_length
Surface roughness in meters for the Global Wind Atlas climate.
The default is the WAsP value for water.
reference_height
Height in meters of the Global Wind Atlas climate. Defaults to the
median hub height of the farm.
Raises
------
KeyError
If ``name`` matches nothing.
TypeError
If ``name``, ``coordinate_frame`` or a non-``None`` ``resource``
is not a string.
ValueError
If ``name`` is ambiguous or only partly matches, or an option is
invalid.
"""
if not isinstance(name, str):
raise TypeError("name must be a string")
if not isinstance(coordinate_frame, str):
raise TypeError("coordinate_frame must be a string")
if resource is not None and not isinstance(resource, str):
raise TypeError("resource must be None or a string")
if resource not in (None, "gwa"):
raise ValueError("resource must be None or 'gwa'")
kind, index = _resolve_name(name)
load = _load_group if kind == "group" else _load_atomic
return load(
index, coordinate_frame, resource, roughness_length, reference_height
)
[docs]
def list(
self,
country: str | None = None,
commissioning_year_range: tuple[int, int] | None = None,
min_turbines: int | None = None,
*,
include_groups: bool = False,
) -> tuple[ReferenceFarmInfo, ...]:
"""List packaged farms without loading their geometry.
``country`` is matched ignoring case, ``commissioning_year_range`` is an
inclusive ``(first, last)`` pair that excludes farms of unknown year,
and ``min_turbines`` is a positive lower bound. Named groups are left
out unless ``include_groups`` is true, so no farm is listed twice.
Raises
------
TypeError
If ``country`` is not a string, ``commissioning_year_range`` is not
a sequence or has an end that is not a number, or ``min_turbines``
is not a number.
ValueError
If ``commissioning_year_range`` does not have two ends, has an end
that is not a whole number or is not ordered, or ``min_turbines``
is not a positive whole number.
"""
if country is not None and not isinstance(country, str):
raise TypeError("country must be a string or None")
if commissioning_year_range is not None:
if isinstance(commissioning_year_range, (str, bytes)) or not hasattr(
commissioning_year_range, "__len__"
):
raise TypeError(
"commissioning_year_range must be a (first, last) pair of years"
)
if len(commissioning_year_range) != 2:
raise ValueError(
"commissioning_year_range must be a (first, last) pair of years"
)
lo, hi = (
_validate_integer("commissioning_year_range", year, minimum=-(2**63))
for year in commissioning_year_range
)
if lo > hi:
raise ValueError("commissioning_year_range must be ordered")
else:
lo = hi = 0
if min_turbines is not None:
min_turbines = _validate_integer("min_turbines", min_turbines)
results = _catalog_infos(include_groups)
return tuple(
replace(
item,
boundary_adjustment=copy.deepcopy(item.boundary_adjustment),
source_dates=dict(item.source_dates),
)
for item in results
if (country is None or item.country.casefold() == country.casefold())
and (min_turbines is None or item.n_turbines >= int(min_turbines))
and (
commissioning_year_range is None
or (
item.commissioning_year is not None
and lo <= item.commissioning_year <= hi
)
)
)
reference_farm = ReferenceCatalog()
