Skip to content

Correct condition 5 in detect_clearksy #596

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 11 commits into from
Oct 9, 2018
34 changes: 23 additions & 11 deletions pvlib/clearsky.py
Original file line number Diff line number Diff line change
Expand Up @@ -680,36 +680,39 @@ 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 +
sample_interval*sample_interval), axis=0)
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 +721,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

Expand Down Expand Up @@ -756,6 +760,14 @@ 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(
Copy link
Member

Choose a reason for hiding this comment

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

maybe "values" instead of "array"?

Copy link
Member Author

Choose a reason for hiding this comment

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

_flag on the previous set of keys, removed _array here

Copy link
Member

Choose a reason for hiding this comment

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

I like the change, but it would break user code if people are calling return_components. Probably not likely to be a serious issue but we should think twice about doing this. In any case it should be listed as an API change.

Copy link
Member Author

Choose a reason for hiding this comment

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

what was meant by "values" instead of "array"? I'm ok reverting this change, but I thought you were requesting renaming the keys.

Copy link
Member

Choose a reason for hiding this comment

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

Sorry for not being clear. I originally meant to suggest only slightly different names for the new keys that you've added. I was not suggesting renaming the existing keys. I thought "array" was an undesirable identifier because everything is an array here.

I do like the latest proposal, but we'd need to document it as an API change.

meas_slope - alpha * clear_slope, axis=0))

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