"""Methods for calculating Shear Exponent (Alpha)"""
import logging
import xarray as xr
import numpy as np
from windkit import spatial as wksp
logger = logging.getLogger(__name__)
ATTRS = {
"long_name": "Wind shear exponent",
"grid_mapping": "crs",
}
[docs]
def shear_exp_from_wspd(wspd, height_map, w_map=None):
"""Calculate shear exponent using a linear log-log fit from a number of heights
.. warning::
This function is experimental and its signature may change.
Parameters
----------
wspd : xr.DataArray
Wind Speeds with at least a height dimension
height_map : dict
Dictionary:
- keys: a float providing the height of Alpha to be calculated
- values: a list of the heights used to calculate Alpha, must be in wspd
w_map : dict
The same as height_map, but values should be the weights that are passed to the
numpy.polyfit function.
Returns
-------
xr.DataArray
DataArray of Alpha
"""
## We don't support 3D point data yet
if wksp.is_point(wspd):
msg = (
"Shear exponent calculation from point objects is not supported yet. "
"If you still want to perform this calculation, convert to stacked_point, "
"or raster format."
)
raise NotImplementedError(msg)
if w_map is not None:
logging.info("Using weights for fitting.")
raise NotImplementedError("Fitting weights are not implemented yet.")
# if w_map.keys != height_map.keys:
# raise ValueError("w_map must have the same keys as height_map.")
# else:
# pass
# Calculate Alpha
hgt_list = []
for hgt, ws_levs in height_map.items():
logger.info("Calculating Alpha at height %s from heights %s", hgt, ws_levs)
# Subset to selected heights
tmp = wspd.sel(height=ws_levs).copy()
# Stack to work with 1D numpy functions
stack_dims = list(tmp.dims)
stack_dims.remove("height") # This modifies list in place returning none
logger.debug("stack_dims: %s", stack_dims)
if len(stack_dims) > 1:
tmp = tmp.stack(tmp_dim=stack_dims).reset_index("tmp_dim")
core_dim = "tmp_dim"
else:
core_dim = stack_dims[0]
fits = xr.apply_ufunc(
np.polyfit,
np.log(tmp.height),
np.log(tmp),
1, # args to polyfit
input_core_dims=[["height"], ["height", core_dim], []],
output_core_dims=[["rank", core_dim]],
)
# Reshape to original format
alpha = fits[0]
if core_dim == "tmp_dim":
alpha = alpha.set_index(tmp_dim=stack_dims).unstack("tmp_dim")
alpha = np.fabs(alpha)
alpha.expand_dims({"height": [hgt]})
alpha.name = "ALPHA"
# Concatenate to list
hgt_list.append(alpha)
logger.info("Concatenating heights to a single DataArray.")
rslt = xr.concat(hgt_list, "height")
rslt.attrs = ATTRS
# Set coords not height to input values
if "stacked_point" in rslt.coords:
del rslt.coords["stacked_point"]
for coord in dict(rslt.coords).copy():
if coord != "height":
rslt.coords[coord] = wspd.coords[coord]
return rslt