Iss494 timeseries maximise data

CHANGELOG.md

@@ -20,7 +20,8 @@ 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))
+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. 

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
@@ -109,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):
@@ -832,27 +841,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
@@ -1094,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]
@@ -1115,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:
@@ -1124,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.
 
@@ -1135,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:
         """
 
@@ -1145,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))
 

tests/test_shear.py

@@ -278,6 +278,100 @@ 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['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
+        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_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']]
@@ -314,15 +408,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)