Skip to content
29 changes: 19 additions & 10 deletions pvlib/clearsky.py
Original file line number Diff line number Diff line change
Expand Up @@ -680,35 +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_ghi_diff = np.diff(measured[H], n=1, axis=0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I prefer meas_diff to meas_ghi_diff. None of the other code in this function specifically refers to ghi.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

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_slope_max = np.max(np.abs(meas_slope), axis=0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

delete this line

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gone

meas_line_length = np.sum(np.sqrt(
meas_slope*meas_slope + sample_interval*sample_interval), axis=0)
meas_ghi_diff*meas_ghi_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_ghi_diff = np.diff(clearsky[H], n=1, axis=0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comment as above re "ghi"

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same action

clear_slope = np.diff(clearsky[H], n=1, axis=0) / sample_interval
# clear_slope_max = np.max(np.abs(clear_slope), axis=0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

delete this line

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gone


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 +
alpha*alpha*clear_ghi_diff*clear_ghi_diff +
sample_interval*sample_interval), axis=0)

line_diff = meas_line_length - clear_line_length
Expand All @@ -718,7 +720,7 @@ 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

Expand Down Expand Up @@ -756,6 +758,13 @@ def rmse(alpha):
components['slope_max'] = c5
components['mean_nan'] = c6
components['windows'] = clear_windows

components['mean_diff_array'] = np.abs(meas_mean - alpha*clear_mean)
components['max_diff_array'] = np.abs(meas_max - alpha*clear_max)
components['line_length_array'] = meas_line_length - clear_line_length
components['slope_nstd_array'] = meas_slope_nstd
components['slope_max_array'] = (np.max(meas_slope - alpha*clear_slope, axis=0))

return clear_samples, components, alpha
else:
return clear_samples
Expand Down