diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst index 7633e2dda1..4c56e91fcd 100644 --- a/docs/sphinx/source/whatsnew/v0.6.1.rst +++ b/docs/sphinx/source/whatsnew/v0.6.1.rst @@ -18,6 +18,7 @@ API Changes and :py:func:`pvlib.iotools.read_tmy3` instead. (:issue:`261`) * Added keyword argument ``horizon`` to :func:`~pvlib.solarposition.pyephem` and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'`` +* Changed key names for `components` returned from :py:func:`pvlib.clearsky.detect_clearsky` Enhancements @@ -35,6 +36,7 @@ Bug fixes ~~~~~~~~~ * Fix when building documentation using Matplotlib 3.0 or greater. * Fix and improve :func:`~pvlib.solarposition.hour_angle` (:issue:`598`) +* Fix error in :func:`pvlib.clearsky.detect_clearsky` (:issue:`506`) Testing @@ -45,6 +47,7 @@ Testing Contributors ~~~~~~~~~~~~ * Will Holmgren (:ghuser:`wholmgren`) -* Cliff Hansen (:ghuser:`cwhanse`) * Leland Boeman (:ghuser:`lboeman`) +* Ben Ellis (:ghuser:`bhellis725`) +* Cliff Hansen (:ghuser:`cwhanse`) * Mark Mikofski (:ghuser:`mikofski`) diff --git a/pvlib/clearsky.py b/pvlib/clearsky.py index 2d85fbe880..be8ca8f46c 100644 --- a/pvlib/clearsky.py +++ b/pvlib/clearsky.py @@ -680,36 +680,37 @@ def detect_clearsky(measured, clearsky, times, window_length, raise NotImplementedError('algorithm does not yet support unequal ' 'times. consider resampling your data.') - samples_per_window = int(window_length / sample_interval) + intervals_per_window = int(window_length / sample_interval) # generate matrix of integers for creating windows with indexing from scipy.linalg import hankel - H = hankel(np.arange(samples_per_window), # noqa: N806 - np.arange(samples_per_window-1, len(times))) + H = hankel(np.arange(intervals_per_window), # noqa: N806 + np.arange(intervals_per_window - 1, len(times))) # calculate measurement statistics meas_mean = np.mean(measured[H], axis=0) meas_max = np.max(measured[H], axis=0) - meas_slope = np.diff(measured[H], n=1, axis=0) + meas_diff = np.diff(measured[H], n=1, axis=0) + meas_slope = np.diff(measured[H], n=1, axis=0) / sample_interval # matlab std function normalizes by N-1, so set ddof=1 here meas_slope_nstd = np.std(meas_slope, axis=0, ddof=1) / meas_mean - meas_slope_max = np.max(np.abs(meas_slope), axis=0) meas_line_length = np.sum(np.sqrt( - meas_slope*meas_slope + sample_interval*sample_interval), axis=0) + meas_diff * meas_diff + + sample_interval * sample_interval), axis=0) # calculate clear sky statistics clear_mean = np.mean(clearsky[H], axis=0) clear_max = np.max(clearsky[H], axis=0) - clear_slope = np.diff(clearsky[H], n=1, axis=0) - clear_slope_max = np.max(np.abs(clear_slope), axis=0) + clear_diff = np.diff(clearsky[H], n=1, axis=0) + clear_slope = np.diff(clearsky[H], n=1, axis=0) / sample_interval from scipy.optimize import minimize_scalar alpha = 1 for iteration in range(max_iterations): clear_line_length = np.sum(np.sqrt( - alpha*alpha*clear_slope*clear_slope + - sample_interval*sample_interval), axis=0) + alpha * alpha * clear_diff * clear_diff + + sample_interval * sample_interval), axis=0) line_diff = meas_line_length - clear_line_length @@ -718,7 +719,8 @@ def detect_clearsky(measured, clearsky, times, window_length, c2 = np.abs(meas_max - alpha*clear_max) < max_diff c3 = (line_diff > lower_line_length) & (line_diff < upper_line_length) c4 = meas_slope_nstd < var_diff - c5 = (meas_slope_max - alpha*clear_slope_max) < slope_dev + c5 = np.max(np.abs(meas_slope - + alpha * clear_slope), axis=0) < slope_dev c6 = (clear_mean != 0) & ~np.isnan(clear_mean) clear_windows = c1 & c2 & c3 & c4 & c5 & c6 @@ -749,13 +751,21 @@ def rmse(alpha): if return_components: components = OrderedDict() - components['mean_diff'] = c1 - components['max_diff'] = c2 - components['line_length'] = c3 - components['slope_nstd'] = c4 - components['slope_max'] = c5 - components['mean_nan'] = c6 + components['mean_diff_flag'] = c1 + components['max_diff_flag'] = c2 + components['line_length_flag'] = c3 + components['slope_nstd_flag'] = c4 + components['slope_max_flag'] = c5 + components['mean_nan_flag'] = c6 components['windows'] = clear_windows + + components['mean_diff'] = np.abs(meas_mean - alpha * clear_mean) + components['max_diff'] = np.abs(meas_max - alpha * clear_max) + components['line_length'] = meas_line_length - clear_line_length + components['slope_nstd'] = meas_slope_nstd + components['slope_max'] = (np.max( + meas_slope - alpha * clear_slope, axis=0)) + return clear_samples, components, alpha else: return clear_samples