Skip to content

Commit

Permalink
Merge pull request #13 from sede-open/optimize_circstd_implementation
Browse files Browse the repository at this point in the history
Optimize wind turbulence calculation
  • Loading branch information
bvandekerkhof authored Sep 12, 2024
2 parents f0c77ef + 8f8ff1e commit 5c1af2b
Showing 3 changed files with 11 additions and 9 deletions.
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
@@ -8,7 +8,7 @@ build-backend = "poetry.core.masonry.api"

[tool.poetry]
name = "pyelq-sdk"
version = "1.0.8b"
version = "1.0.9"
description = "Package for detection, localization and quantification code."
authors = ["Bas van de Kerkhof", "Matthew Jones", "David Randell"]
homepage = "https://sede-open.github.io/pyELQ/"
14 changes: 9 additions & 5 deletions src/pyelq/meteorology.py
Original file line number Diff line number Diff line change
@@ -16,7 +16,6 @@
import plotly.express as px
import plotly.graph_objects as go
from pandas.arrays import DatetimeArray
from scipy.stats import circstd

from pyelq.coordinate_system import Coordinate
from pyelq.sensor.sensor import SensorGroup
@@ -102,7 +101,12 @@ def calculate_uv_from_wind_speed_direction(self) -> None:
def calculate_wind_turbulence_horizontal(self, window: str) -> None:
"""Calculate the horizontal wind turbulence values from the wind direction attribute.
Wind turbulence values are calculated as the circular standard deviation based on a rolling window.
Wind turbulence values are calculated as the circular standard deviation of wind direction
(https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.circstd.html).
The implementation here is equivalent to using the circstd function from scipy.stats as an apply
function on a rolling window. However, using the rolling mean on sin and cos speeds up
the calculation by a factor of 100.
Outputted values are calculated at the center of the window and at least 3 observations are required in a
window for the calculation. If the window contains less values the result will be np.nan.
The result of the calculation will be stored as the wind_turbulence_horizontal attribute.
@@ -113,9 +117,9 @@ def calculate_wind_turbulence_horizontal(self, window: str) -> None:
"""
data_series = pd.Series(data=self.wind_direction, index=self.time)
aggregated_data = data_series.rolling(window=window, center=True, min_periods=3).apply(
circstd, kwargs={"low": 0, "high": 360}
)
sin_rolling = (np.sin(data_series * np.pi / 180)).rolling(window=window, center=True, min_periods=3).mean()
cos_rolling = (np.cos(data_series * np.pi / 180)).rolling(window=window, center=True, min_periods=3).mean()
aggregated_data = np.sqrt(-2 * np.log((sin_rolling**2 + cos_rolling**2) ** 0.5)) * 180 / np.pi
self.wind_turbulence_horizontal = aggregated_data.values

def plot_polar_hist(self, nof_sectors: int = 16, nof_divisions: int = 5, template: object = None) -> go.Figure():
4 changes: 1 addition & 3 deletions tests/test_meteorology.py
Original file line number Diff line number Diff line change
@@ -201,7 +201,7 @@ def test_meteorology_group():
def test_calculate_wind_turbulence_horizontal():
"""Checks that the wind turbulence values are calculated correctly.
To verify circstd, we define winds as draws from a normal distribution. We then check that the mean of the
To verify horizontal wind turbulence calculations, we define winds as draws from a normal distribution. We then check that the mean of the
calculated turbulence values is within 3 standard deviations of the true value.
"""
@@ -219,10 +219,8 @@ def test_calculate_wind_turbulence_horizontal():
pd.date_range(dt.datetime(2023, 1, 1), dt.datetime(2023, 1, 2), freq="5s"), dtype="datetime64[ns]"
)
met.wind_direction = np.random.normal(180, sigma, met.time.shape[0])

met.calculate_wind_turbulence_horizontal(window="300s")

tolerance = 3 * np.std(met.wind_turbulence_horizontal)
mean_turbulence = np.mean(met.wind_turbulence_horizontal)

assert (mean_turbulence - tolerance) < sigma < (mean_turbulence + tolerance)

0 comments on commit 5c1af2b

Please sign in to comment.