"""Align wind-climate direction frames between coordinate reference systems."""
import numpy as np
import pyproj
import xarray as xr
import windkit as wk
from pywasp._batch import batch_layout
from pywasp._windkit_internal import _spatial_dims_in_order
from pywasp.wasp.binned_wind_climate import bwc_resample_sectors
from pywasp.wasp.time_series_wind_climate import tswc_rotate
from pywasp.wasp.weibull_wind_climate import wwc_rotate
__all__ = ["align_direction_crs"]
_WIND_DIR_CRS_CAPABILITIES = (
("geowc", wk.is_geowc),
("tswc", wk.is_tswc),
("bwc", wk.is_bwc),
("wwc", wk.is_wwc),
("gwc", wk.is_gwc),
)
[docs]
def align_direction_crs(
wind_climate,
target_wind_dir_crs,
source_wind_dir_crs=None,
conf=None,
):
"""Align a wind climate to a target wind-direction CRS.
Wind directions may be expressed relative to true north or relative to the
grid north of a projected CRS. This helper rotates the directional part of a
wind climate from ``source_wind_dir_crs`` into ``target_wind_dir_crs`` at the
climate locations and stores the target frame in the ``wind_dir_crs``
attribute.
The simpler the wind climate, the harder this operation is to perform
without information loss. TSWC direction values are rotated directly and
therefore retain the original time-series information. BWC sector histograms
are resampled, so information is limited by sector and wind-speed bin
resolution. WWC and GWC contain only sectoral Weibull parameters and
frequencies, so rotation is the most compressed approximation and relies on
conserving moments during sector rotation.
Parameters
----------
wind_climate : xarray.Dataset
TSWC, BWC, WWC, or GWC dataset. GeoWC is not supported because it
contains additional directional variables that cannot be rotated by the
WWC sector-rotation routine.
target_wind_dir_crs : CRS-like
Coordinate reference system whose grid north should define the output
wind-direction frame.
source_wind_dir_crs : CRS-like, optional
Coordinate reference system whose grid north defines the input
wind-direction frame. If omitted, the function reads
``wind_climate.attrs["wind_dir_crs"]``. All supported wind climates
require explicit source-frame metadata.
conf : pywasp.wasp.Config, optional
Configuration object passed to ``wwc_rotate`` for WWC and GWC input.
Returns
-------
xarray.Dataset
Wind climate with directions aligned to ``target_wind_dir_crs`` and a
``wind_dir_crs`` attribute containing the target CRS WKT. The spatial
CRS metadata is preserved from the input.
Raises
------
ValueError
If ``wind_climate`` is unsupported, if source-frame metadata is missing,
or if a GeoWC is passed.
"""
target_crs = pyproj.CRS.from_user_input(target_wind_dir_crs)
_wind_dir_crs_capability(wind_climate)
source_crs = _source_wind_dir_crs(
wind_climate,
source_wind_dir_crs=source_wind_dir_crs,
)
if source_crs.equals(target_crs):
result = wind_climate.copy()
else:
rotation = _direction_frame_rotation(wind_climate, target_crs, source_crs)
result = _rotate_wind_climate(
wind_climate,
rotation_angle_deg=-rotation,
conf=conf,
)
result.attrs["wind_dir_crs"] = target_crs.to_wkt()
return result
def _source_wind_dir_crs(wind_climate, source_wind_dir_crs):
"""Return the source wind-direction CRS from input or metadata.
Parameters
----------
wind_climate : xarray.Dataset
Wind climate dataset that may contain a ``wind_dir_crs`` attribute.
source_wind_dir_crs : CRS-like or None
Explicit source wind-direction CRS.
Returns
-------
pyproj.CRS
Source wind-direction CRS.
Raises
------
ValueError
If no source CRS is provided and ``wind_climate`` has no
``wind_dir_crs`` attribute.
"""
if source_wind_dir_crs is not None:
return pyproj.CRS.from_user_input(source_wind_dir_crs)
wind_dir_crs = wind_climate.attrs.get("wind_dir_crs")
if wind_dir_crs is not None:
return pyproj.CRS.from_user_input(wind_dir_crs)
raise ValueError(
"source_wind_dir_crs is required when wind_climate has no wind_dir_crs "
"attribute."
)
def _direction_frame_rotation(wind_climate, target_crs, source_crs):
"""Return target-frame grid convergence on the input spatial structure.
Grid convergence only depends on horizontal position, so the returned
rotation carries no independent ``height`` dimension. The rotation
consumers do not yet broadcast it over multi-height cuboid or stacked-point
climates; that pre-existing limitation is tracked in #757. Restoration
uses the input's exact coordinates because a CRS round trip introduces
floating-point coordinate noise. Internal reprojection does not replace the
input spatial CRS metadata.
Parameters
----------
wind_climate : xarray.Dataset
Wind climate dataset whose spatial structure should be preserved.
target_crs : pyproj.CRS
Target wind-direction CRS.
source_crs : pyproj.CRS
Source wind-direction CRS.
Returns
-------
xarray.DataArray
Grid-convergence rotation on the original wind-climate spatial
structure.
"""
horizontal = wind_climate[["west_east", "south_north"]]
layout = batch_layout(
horizontal,
batch_dims=_spatial_dims_in_order(horizontal),
core_dims=(),
required_coords=("west_east", "south_north"),
)
source_we, source_sn = layout.flat_coords(dtype=np.float64)
# Use WindKit's direct reprojection once windkit#811 preserves the input
# raster/cuboid structure instead of restacking it onto a point dimension.
transformer = pyproj.Transformer.from_crs(
wk.spatial.get_crs(wind_climate), target_crs, always_xy=True
)
target_we, target_sn = transformer.transform(source_we, source_sn)
target_points = wk.spatial.set_crs(
xr.Dataset(
{"west_east": ("_flat", target_we), "south_north": ("_flat", target_sn)}
).set_coords(["west_east", "south_north"]),
target_crs,
)
rotation = wk.spatial.grid_convergence(target_points, source=source_crs)
return layout.restore(
rotation.data, core_dims=(), name=rotation.name, attrs=rotation.attrs
)
def _rotate_wind_climate(
wind_climate,
rotation_angle_deg,
conf,
):
"""Rotate a wind climate using the climate-specific implementation.
Parameters
----------
wind_climate : xarray.Dataset
TSWC, BWC, WWC, or GWC dataset.
rotation_angle_deg : float or xarray.DataArray
Clockwise-positive rotation angle in degrees.
conf : pywasp.wasp.Config or None
Configuration object passed to ``wwc_rotate`` for WWC and GWC input.
Returns
-------
xarray.Dataset
Rotated wind climate dataset.
"""
capability = _wind_dir_crs_capability(wind_climate)
if capability == "tswc":
return tswc_rotate(wind_climate, rotation_angle_deg)
if capability == "bwc":
return bwc_resample_sectors(
wind_climate,
n_sectors=wind_climate.sizes["sector"],
offset=rotation_angle_deg,
)
if capability in ("wwc", "gwc"):
return wwc_rotate(wind_climate, rotation_angle_deg, conf=conf)
raise ValueError(
"align_direction_crs only supports TSWC, BWC, WWC, and GWC datasets."
)
def _wind_dir_crs_capability(wind_climate):
"""Return the supported wind-climate capability for direction alignment.
Parameters
----------
wind_climate : xarray.Dataset
Wind climate dataset to classify.
Returns
-------
str
One of ``"tswc"``, ``"bwc"``, ``"wwc"``, or ``"gwc"``.
Raises
------
ValueError
If ``wind_climate`` is GeoWC or unsupported.
"""
for name, check in _WIND_DIR_CRS_CAPABILITIES:
if check(wind_climate):
if name == "geowc":
raise ValueError(
"GeoWC is not supported by align_direction_crs because its "
"geostrophic variables cannot be rotated by the WWC sector "
"rotation routine."
)
return name
raise ValueError(
"align_direction_crs only supports TSWC, BWC, WWC, and GWC datasets."
)
