From c8f14ac0473e4c356aea7485dd09f7d3bc23ca6e Mon Sep 17 00:00:00 2001 From: Olivia Bentley Date: Wed, 30 Sep 2026 17:43:38 +0100 Subject: [PATCH 1/2] iss #494 bug fix for maximise data --- CHANGELOG.md | 2 +- brightwind/analyse/shear.py | 56 +++++++++++++++++++-------- tests/test_shear.py | 76 +++++++++++++++++++++++++++++++++++-- 3 files changed, 115 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e309578d..ba4839b6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,7 +20,7 @@ XX-Xxx-2026 1. `Shear.TimeSeries` and `Shear.BySector` no longer print "This may take a while..." and no longer clear the notebook cell output when they run. ([#639](https://github.com/brightwind-dev/brightwind/issues/639)) ### Bug Fixes -1. +1. Fixed `maximise_data=True` having no effect in `Shear.TimeSeries`. Shear is now calculated for every timestamp with at least two valid wind speeds (above `min_speed` and not NaN), as described in the docstring, using only the valid heights at each timestamp. Results with `maximise_data=False` (the default), and for timestamps where all heights are valid, are unchanged. ([#494](https://github.com/brightwind-dev/brightwind/issues/494)) ### Deprecated 1. diff --git a/brightwind/analyse/shear.py b/brightwind/analyse/shear.py index 17f5801b..5f181c96 100644 --- a/brightwind/analyse/shear.py +++ b/brightwind/analyse/shear.py @@ -85,19 +85,27 @@ def __init__(self, wspds, heights, min_speed=3, calc_method='power_law', max_plo wspds, cvg = Shear._data_prep( wspds=wspds, heights=heights, min_speed=min_speed, maximise_data=maximise_data, return_raw_wspds=True ) - valid_mask = (wspds > min_speed).all(axis=1) & ~wspds.isna().any(axis=1) + valid_heights = (wspds > min_speed) & wspds.notna() + if maximise_data: + # use every timestamp with at least two valid heights, fitting to only the valid heights + valid_mask = valid_heights.sum(axis=1) >= 2 + fit_mask = valid_heights.loc[valid_mask].values + else: + valid_mask = valid_heights.all(axis=1) + fit_mask = None + valid_wspds = wspds.loc[valid_mask].where(valid_heights.loc[valid_mask]).values log_heights = np.log(np.asarray(heights, dtype=float)) if calc_method == 'power_law': alpha = pd.Series(index=wspds.index, dtype=float, name='alpha') - alpha.loc[valid_mask], _ = Shear._calc_linear_fit(log_heights, np.log(wspds.loc[valid_mask].values)) + alpha.loc[valid_mask], _ = Shear._calc_linear_fit(log_heights, np.log(valid_wspds), mask=fit_mask) self._alpha = alpha elif calc_method == 'log_law': slope = pd.Series(index=wspds.index, dtype=float) intercept = pd.Series(index=wspds.index, dtype=float) slope.loc[valid_mask], intercept.loc[valid_mask] = Shear._calc_linear_fit( - log_heights, wspds.loc[valid_mask].values) + log_heights, valid_wspds, mask=fit_mask) roughness_coefficient = pd.Series(Shear._calc_roughness(slope=slope, intercept=intercept), name='roughness_coefficient') self._roughness = roughness_coefficient @@ -832,27 +840,45 @@ def _calc_power_law(wspds, heights, return_coeff=False, maximise_data=False): return coeffs[0] @staticmethod - def _calc_linear_fit(x, y): + def _calc_linear_fit(x, y, mask=None): """ Derive the least squares straight line fit of y against x for every row of y at once. This gives the same result as calling np.polyfit(x, y_row, deg=1) on each row, but is vectorised across all rows. - :param x: Values to fit against, e.g. log of heights, one per column of y. - :type x: numpy.ndarray - :param y: 2D array with one row per timestamp and one column per height, e.g. log of wind speeds. Each row - must be complete (no NaN). - :type y: numpy.ndarray - :return: The slope and intercept of the best fit line for each row. - :rtype: tuple(numpy.ndarray, numpy.ndarray) + :param x: Values to fit against, e.g. log of heights, one per column of y. + :type x: numpy.ndarray + :param y: 2D array with one row per timestamp and one column per height, e.g. log of wind speeds. If mask + is None each row must be complete (no NaN). + :type y: numpy.ndarray + :param mask: Optional 2D boolean array, the same shape as y, which is True for the values of y to use in the + fit of each row. Rows need at least two values with different x, otherwise the slope and + intercept are NaN. + :type mask: numpy.ndarray or None + :return: The slope and intercept of the best fit line for each row. + :rtype: tuple(numpy.ndarray, numpy.ndarray) METHODOLOGY: slope = sum((x - mean(x)) * (y - mean(y))) / sum((x - mean(x))^2) intercept = mean(y) - slope * mean(x) + where the sums and means of each row only use the values selected by mask. """ - x_dev = x - x.mean() - y_mean = y.mean(axis=1) - slope = ((y - y_mean[:, np.newaxis]) * x_dev).sum(axis=1) / (x_dev ** 2).sum() - intercept = y_mean - slope * x.mean() + if mask is None: + x_dev = x - x.mean() + y_mean = y.mean(axis=1) + slope = ((y - y_mean[:, np.newaxis]) * x_dev).sum(axis=1) / (x_dev ** 2).sum() + intercept = y_mean - slope * x.mean() + return slope, intercept + + x = np.broadcast_to(x, y.shape) + count = mask.sum(axis=1) + with np.errstate(divide='ignore', invalid='ignore'): + x_mean = np.where(mask, x, 0).sum(axis=1) / count + y_mean = np.where(mask, y, 0).sum(axis=1) / count + x_dev = np.where(mask, x - x_mean[:, np.newaxis], 0) + y_dev = np.where(mask, y - y_mean[:, np.newaxis], 0) + x_dev_sq_sum = (x_dev ** 2).sum(axis=1) + slope = np.where(x_dev_sq_sum > 0, (x_dev * y_dev).sum(axis=1) / x_dev_sq_sum, np.nan) + intercept = y_mean - slope * x_mean return slope, intercept @staticmethod diff --git a/tests/test_shear.py b/tests/test_shear.py index 680d19f9..f6722691 100644 --- a/tests/test_shear.py +++ b/tests/test_shear.py @@ -278,6 +278,70 @@ def test_calc_linear_fit(): np.testing.assert_allclose(slope, expected_slope, rtol=0, atol=1e-12) np.testing.assert_allclose(intercept, expected_intercept, rtol=0, atol=1e-12) + # fit only the values selected by the mask in each row + heights = [40, 60, 80, 100, 120] + log_heights = np.log(heights) + log_wspds = np.log(rng.uniform(3, 20, size=(200, len(heights)))) + mask = rng.random(size=log_wspds.shape) > 0.4 + mask[:, :2] = [True, True] + log_wspds[~mask] = np.nan + slope, intercept = bw.Shear._calc_linear_fit(log_heights, log_wspds, mask=mask) + for i in range(len(log_wspds)): + expected_slope, expected_intercept = np.polyfit(log_heights[mask[i]], log_wspds[i, mask[i]], deg=1) + assert slope[i] == pytest.approx(expected_slope, abs=1e-12) + assert intercept[i] == pytest.approx(expected_intercept, abs=1e-12) + # a full mask gives the same result as no mask + full_slope, full_intercept = bw.Shear._calc_linear_fit(log_heights, log_wspds[mask.all(axis=1)]) + np.testing.assert_array_equal(slope[mask.all(axis=1)], full_slope) + # a row with fewer than two different heights can't be fitted + slope, intercept = bw.Shear._calc_linear_fit(np.log([80, 80, 60]), np.log([[8.0, 8.1, 7.5], [8.0, np.nan, 7.5]]), + mask=np.array([[True, True, False], [True, False, False]])) + assert np.isnan(slope).all() and np.isnan(intercept).all() + + +def test_time_series_maximise_data(): + heights = [80, 60, 40] + # NaN in one height and below min_speed in another, so some timestamps only have two valid heights + anemometers = DATA[['Spd80mN', 'Spd60mN', 'Spd40mN']].copy() + anemometers.iloc[::5, 0] = np.nan + anemometers.iloc[1::7, 1] = 2.0 + anemometers.iloc[2::11, 2] = np.nan + valid = (anemometers > 3).values + two_valid = valid.sum(axis=1) == 2 + # positions of the two valid heights in each row with exactly two + first, second = np.argsort(~valid[two_valid], axis=1, kind='stable')[:, :2].T + speeds = anemometers.values[two_valid] + rows = np.arange(len(speeds)) + u1, u2 = speeds[rows, first], speeds[rows, second] + z1, z2 = np.array(heights)[first], np.array(heights)[second] + + for calc_method in ['power_law', 'log_law']: + all_heights = bw.Shear.TimeSeries(anemometers, heights, calc_method=calc_method) + maximised = bw.Shear.TimeSeries(anemometers, heights, calc_method=calc_method, maximise_data=True) + name = 'alpha' if calc_method == 'power_law' else 'roughness' + shear_all, shear_max = getattr(all_heights, name), getattr(maximised, name) + + assert shear_max.index.equals(anemometers.index) + assert all_heights.info['output data']['concurrent_period'] == 0.941 + assert maximised.info['output data']['concurrent_period'] == 1.455 + assert shear_all.notna().sum() == 49513 + assert shear_max.notna().sum() == 76552 == (valid.sum(axis=1) >= 2).sum() + # timestamps with all heights valid are unchanged, timestamps with fewer than two are NaN + complete = valid.all(axis=1) + np.testing.assert_array_equal(shear_max[complete].values, shear_all[complete].values) + assert shear_max[valid.sum(axis=1) < 2].isna().all() + + # timestamps with two valid heights use only those two + if calc_method == 'power_law': + expected = np.log(u1 / u2) / np.log(z1 / z2) + np.testing.assert_allclose(shear_max[two_valid].values, expected, rtol=0, atol=1e-12) + else: + slope = (u1 - u2) / np.log(z1 / z2) + intercept = u1 - slope * np.log(z1) + expected = np.exp(-intercept / slope) + realistic = expected < 10 + np.testing.assert_allclose(shear_max[two_valid].values[realistic], expected[realistic], rtol=1e-8) + def test_time_series_full_data(): anemometers = DATA[['Spd80mN', 'Spd60mN', 'Spd40mN']] @@ -314,15 +378,21 @@ def test_time_series(): anemometers = anemometers[:100] # Test initialisation shear_by_ts_power_law = bw.Shear.TimeSeries(anemometers, heights) - shear_by_ts_power_law = bw.Shear.TimeSeries(anemometers, heights, maximise_data=True) shear_by_ts_log_law = bw.Shear.TimeSeries(anemometers, heights, calc_method='log_law') - shear_by_ts_log_law = bw.Shear.TimeSeries(anemometers, heights, calc_method='log_law', - maximise_data=True) + shear_by_ts_power_law_max = bw.Shear.TimeSeries(anemometers, heights, maximise_data=True) + shear_by_ts_log_law_max = bw.Shear.TimeSeries(anemometers, heights, calc_method='log_law', + maximise_data=True) # Test attributes assert round(shear_by_ts_power_law.alpha.mean(), 4) == 0.1786 # Changed to support equality for very large numbers assert abs(shear_by_ts_log_law.roughness.mean() / 4.306534305567819e+68 - 1) < 1e-6 + # maximise_data also uses timestamps where only two heights are valid + assert shear_by_ts_power_law.alpha.count() == 90 + assert shear_by_ts_power_law_max.alpha.count() == 93 + assert shear_by_ts_power_law_max.alpha.mean() == pytest.approx(0.191743, abs=1e-6) + assert shear_by_ts_log_law_max.roughness.count() == 93 + assert shear_by_ts_log_law_max.roughness.median() == pytest.approx(0.562503, abs=1e-6) # Test apply shear_by_ts_power_law.apply(DATA['Spd80mN'], 40, 60) From 8c37d6f7b79b22d990477b7295548619a8daa8c8 Mon Sep 17 00:00:00 2001 From: Olivia Bentley Date: Thu, 22 Oct 2026 16:10:54 +0100 Subject: [PATCH 2/2] iss #494 fix maximise_data issue --- CHANGELOG.md | 1 + brightwind/analyse/shear.py | 20 ++++++++++++++++---- tests/test_shear.py | 30 ++++++++++++++++++++++++++++++ 3 files changed, 47 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ba4839b6..6e5bad40 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ XX-Xxx-2026 ### Bug Fixes 1. Fixed `maximise_data=True` having no effect in `Shear.TimeSeries`. Shear is now calculated for every timestamp with at least two valid wind speeds (above `min_speed` and not NaN), as described in the docstring, using only the valid heights at each timestamp. Results with `maximise_data=False` (the default), and for timestamps where all heights are valid, are unchanged. ([#494](https://github.com/brightwind-dev/brightwind/issues/494)) +1. `Shear.TimeSeries` with `maximise_data=True` no longer fails its data check when wind speeds are never valid at all heights at the same timestamp, or when one height has no valid data. It now only needs at least one timestamp with two or more valid wind speeds. `maximise_data` is now also recorded in `Shear.TimeSeries.info`. The wind profile plot still only uses timestamps where all heights are valid, so the full `Shear.TimeSeries` call can still fail on such data until [#341](https://github.com/brightwind-dev/brightwind/issues/341) is addressed. ([#494](https://github.com/brightwind-dev/brightwind/issues/494)) ### Deprecated 1. diff --git a/brightwind/analyse/shear.py b/brightwind/analyse/shear.py index 5f181c96..e77186f5 100644 --- a/brightwind/analyse/shear.py +++ b/brightwind/analyse/shear.py @@ -117,7 +117,8 @@ def __init__(self, wspds, heights, min_speed=3, calc_method='power_law', max_plo self.calc_method = calc_method self.wspds = wspds self.plot = avg_plot.plot - self.info = Shear._create_info(self, heights=heights, cvg=cvg, min_speed=min_speed) + self.info = Shear._create_info(self, heights=heights, cvg=cvg, min_speed=min_speed, + maximise_data=maximise_data) @property def alpha(self): @@ -1120,7 +1121,7 @@ def _data_prep(wspds, heights, min_speed, maximise_data=False, return_raw_wspds= Shear._unequal_wspd_heights_error_msg(wspds, heights) cvg = None - Shear._valid_wsp_data_error_msg(wspds, min_speed) + Shear._valid_wsp_data_error_msg(wspds, min_speed, maximise_data=maximise_data) if isinstance(wspds.index, pd.DatetimeIndex): if maximise_data is False: cvg = coverage(wspds[wspds > min_speed].dropna(), period='1YS').sum().iloc[1] @@ -1141,7 +1142,15 @@ def _unequal_wspd_heights_error_msg(wspds, heights): str(len(heights)) + ' height(s) were given.') @staticmethod - def _valid_wsp_data_error_msg(wspds, min_speed): + def _valid_wsp_data_error_msg(wspds, min_speed, maximise_data=False): + if maximise_data: + # only need at least two valid heights at some timestamp + if (wspds.notna().sum(axis=1) >= 2).sum() == 0: + raise ValueError("There is no valid data within the dataset provided to calculate the shear.") + if ((wspds > min_speed).sum(axis=1) >= 2).sum() == 0: + raise ValueError(f"There is no valid data above {min_speed} m/s within the dataset provided to " + f"calculate the shear.") + return if wspds.isnull().all().any(): raise ValueError("There is no valid data within the dataset provided to calculate the shear.") if len(wspds[wspds > min_speed].dropna()) == 0: @@ -1150,7 +1159,7 @@ def _valid_wsp_data_error_msg(wspds, min_speed): @staticmethod def _create_info(self, heights, min_speed, cvg, direction_bin_array=None, segments_per_day=None, - segment_start_time=None): + segment_start_time=None, maximise_data=None): """ Create the info dict that is returned to the user. @@ -1161,6 +1170,7 @@ def _create_info(self, heights, min_speed, cvg, direction_bin_array=None, segmen :param direction_bin_array: :param segments_per_day: :param segment_start_time: + :param maximise_data: :return: """ @@ -1171,6 +1181,8 @@ def _create_info(self, heights, min_speed, cvg, direction_bin_array=None, segmen 'min_spd': min_speed} input_data['input_wind_speeds'] = input_wind_speeds input_data['calculation_method'] = self.calc_method + if maximise_data is not None: + input_data['maximise_data'] = maximise_data if cvg is not None: output_data['concurrent_period'] = float("{:.3f}".format(cvg)) diff --git a/tests/test_shear.py b/tests/test_shear.py index f6722691..511d3a4f 100644 --- a/tests/test_shear.py +++ b/tests/test_shear.py @@ -322,6 +322,8 @@ def test_time_series_maximise_data(): shear_all, shear_max = getattr(all_heights, name), getattr(maximised, name) assert shear_max.index.equals(anemometers.index) + assert all_heights.info['input data']['maximise_data'] is False + assert maximised.info['input data']['maximise_data'] is True assert all_heights.info['output data']['concurrent_period'] == 0.941 assert maximised.info['output data']['concurrent_period'] == 1.455 assert shear_all.notna().sum() == 49513 @@ -343,6 +345,34 @@ def test_time_series_maximise_data(): np.testing.assert_allclose(shear_max[two_valid].values[realistic], expected[realistic], rtol=1e-8) +def test_time_series_maximise_data_validation(): + heights = [80, 60, 40] + # heights never all valid at the same timestamp + anemometers = DATA[['Spd80mN', 'Spd60mN', 'Spd40mN']].iloc[:100].copy() + anemometers.iloc[::2, 0] = np.nan + anemometers.iloc[1::2, 2] = np.nan + with pytest.raises(ValueError) as except_info: + bw.Shear.TimeSeries(anemometers, heights) + assert str(except_info.value) == "There is no valid data above 3 m/s within the dataset provided to calculate the shear." + + # a whole height NaN is only allowed with maximise_data + with pytest.raises(ValueError) as except_info: + bw.Shear.TimeSeries(anemometers.assign(Spd40mN=np.nan), heights) + assert str(except_info.value) == "There is no valid data within the dataset provided to calculate the shear." + + # maximise_data still needs two valid heights at some timestamp + one_height = anemometers.assign(Spd60mN=np.nan, Spd40mN=np.nan) + with pytest.raises(ValueError) as except_info: + bw.Shear._valid_wsp_data_error_msg(one_height, 3, maximise_data=True) + assert str(except_info.value) == "There is no valid data within the dataset provided to calculate the shear." + below_min_speed = anemometers.assign(Spd60mN=2.0) + with pytest.raises(ValueError) as except_info: + bw.Shear._valid_wsp_data_error_msg(below_min_speed, 3, maximise_data=True) + assert str(except_info.value) == "There is no valid data above 3 m/s within the dataset provided to calculate the shear." + bw.Shear._valid_wsp_data_error_msg(anemometers, 3, maximise_data=True) + bw.Shear._valid_wsp_data_error_msg(anemometers.assign(Spd40mN=np.nan), 3, maximise_data=True) + + def test_time_series_full_data(): anemometers = DATA[['Spd80mN', 'Spd60mN', 'Spd40mN']] heights = [80, 60, 40]