Townsend snow #1251
Townsend snow #1251
Conversation
There was a problem hiding this comment.
Thanks @abhisheksparikh this is a great start. Some comments for you to consider as you work through the formatting.
pvlib/snow.py
Outdated
| ''' | ||
| return(np.where(N>0, 0.5 * S * (1 + 1/N), 0)) | ||
|
|
||
| def townsend_snow_loss_model(S, N, tilt, RH, T_air, POA, R, H, P=40): |
There was a problem hiding this comment.
Prefer loss_townsend, and would prefer more descriptive parameter names e.g., snow_total instead of S, num_events instead of N. It's a good idea to maintain traceability to the terms used in the reference but we can do that in the description of each parameter.
|
Excuse me for the long delay! I will be working on resolving the issues of this PR in the coming couple of days. |
There was a problem hiding this comment.
Thanks @abhisheksparikh, looking good so far! Can you please also add an entry to docs/sphinx/source/whatsnew/v0.9.1.rst under Enhancements, and yourself to the contributors list?
Several other small changes - variable names change, comment change - in response to Kevin's review notes
| @@ -28,6 +28,8 @@ Enhancements | |||
| * Added ``map_variables`` option to :func:`~pvlib.iotools.read_crn` (:pull:`1368`) | |||
| * Added :py:func:`pvlib.temperature.prilliman` for modeling cell temperature | |||
| at short time steps (:issue:`1081`, :pull:`1391`) | |||
| * Added Townsend Powers Snow loss model in :py:func:`pvlib.snow` | |||
There was a problem hiding this comment.
| * Added Townsend Powers Snow loss model in :py:func:`pvlib.snow` | |
| * Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.townsend` |
| ------- | ||
| loss : numeric | ||
| Average monthly DC capacity loss due to snow coverage [%] | ||
|
|
There was a problem hiding this comment.
Let's add a Notes section here pointing out that this model has not been validated for tracking arrays, but the reference suggests using the maximum rotation angle in place of surface_tilt.
| lower_edge_drop_height : float | ||
| Drop height from array edge to ground [in] | ||
|
|
||
| P : float |
There was a problem hiding this comment.
This P slipped through the cracks :)
|
|
||
| GIT = 1 - C2 * np.exp(-gamma) | ||
| loss = (C1 * Se_weighted * (cosd(surface_tilt))**2 * GIT * | ||
| relative_humidity / (temp_air+273.15)**2 / poa_global**0.67) / 100 |
There was a problem hiding this comment.
Good catch that temperature is in Kelvin 👍
| Relative humidity [%] | ||
|
|
||
| temp_air : numeric | ||
| Ambient temperature [°C] |
There was a problem hiding this comment.
I think it's worth explicitly using the word "monthly" for all the monthly inputs, and "average/total/etc" as appropriate. Temperature is monthly average, and I'm not seeing it in the reference but I guess RH is monthly average as well?
| Row length in the slanted plane of array dimension [in] | ||
|
|
||
| lower_edge_drop_height : float | ||
| Drop height from array edge to ground [in] |
There was a problem hiding this comment.
To discuss: what units should snow_total, slant_height, and lower_edge_drop_height use? I propose cm, m, m respectively, all for consistency with other pvlib functions.
There was a problem hiding this comment.
I prefer metric system over everything else :)
There was a problem hiding this comment.
I can change slant_height and lower_edge_drop_height inputs to m. However, I think, we should leave snow_total in in because that's what is available in most of the weather sources. Thoughts?
There was a problem hiding this comment.
Other than a general reluctance to introduce imperial units to pvlib, the main reason I suggested cm is because snow.fully_covered_nrel takes snowfall in cm. Do non-US datasets use inches? I've only ever used snowfall data from the GHCN, which I think reports in inches.
There was a problem hiding this comment.
I'm OK asking users to do some unit conversions to keep consistent units for similar functions in pvlib
There was a problem hiding this comment.
NOAA NCEI provides monthly climate normals in both metric and imperial units.
I agree that sticking with metric units in pvlib is preferable. I think sticking with a single unit system within a single function is a must.
| Returns | ||
| ------- | ||
| loss : numeric | ||
| Average monthly DC capacity loss due to snow coverage [%] |
There was a problem hiding this comment.
| Average monthly DC capacity loss due to snow coverage [%] | |
| Average monthly DC capacity loss fraction due to snow coverage |
| Se_weighted = 1/3 * Se_prev + 2/3 * Se | ||
| gamma = (slant_height * Se_weighted * cosd(surface_tilt)) / \ | ||
| (np.clip((lower_edge_drop_height**2 - Se_weighted**2), a_min=0.01, | ||
| a_max=None) / 2 / tand(angle_of_repose)) |
There was a problem hiding this comment.
I don't disagree with the spirit of having a lower bound of 0.01, but is there anything special about that value other than being small but nonzero? May be worth a comment explaining it.
There was a problem hiding this comment.
I think a_min is set at 0.01 to prevent gamma from blowing up.
There was a problem hiding this comment.
Right, I was more wondering why specifically 0.01 is the lower limit (why not 0.1, 0.001, or any other small nonzero value?). I guess I was more curious than anything else -- the specific lower bound, so long as it's <<1, doesn't seem to matter much. Feel free to ignore this comment :)
There was a problem hiding this comment.
A common small value is eps or sometimes EPS**(1/3) but I don’t know the mathematics behind why. EPS is the smallest number that can represent on the runtime machine. See np.finfo
https://numpy.org/doc/stable/reference/generated/numpy.finfo.html
There was a problem hiding this comment.
In this case, it seems the danger is having negative gamma, which upstream could result in a performance gain (from negative GIT in Eq. 3). To me, clipping gamma to be positive (not zero) is reasonable. It is not clear if gamma should have an upper bound, I haven't untangled the model equations enough to determine,
There was a problem hiding this comment.
This line is pretty difficult for me to parse - can we break it up in a couple of assignments or be more aggressive with line breaks?
There was a problem hiding this comment.
In this case, it seems the danger is having negative gamma, which upstream could result in a performance gain (from negative GIT in Eq. 3). To me, clipping gamma to be positive (not zero) is reasonable.
If I understand correctly... I mostly agree, but why can't gamma be 0? Shouldn't that be expected when snow fall is 0? I think mathematically gamma will blow up when both the snow fall is 0 and the drop height goes to 0. But I think in practice python will evaluate 0 / 0 and return nan. Assuming we want to avoid that... Drop height should always be greater than 0, even if just by 1 mm, so maybe that's the point to insert our judgement into the code.
There was a problem hiding this comment.
My comment above is not correct. Sorry. Gamma will blow up as effective snow fall approaches the drop height. Beyond that, gamma is large and negative, eventually approaching 0 again.
The ground interference term goes negative when gamma is more negative than ln(c2) ~= -0.67.
|
|
||
| Parameters | ||
| ---------- | ||
| snow_total : numeric |
There was a problem hiding this comment.
Our current definition of numeric includes scalars (e.g. float), which I don't think makes sense here. Would array-like make sense? See https://pvlib-python.readthedocs.io/en/stable/contributing.html#documentation
| np.testing.assert_allclose(expected, actual, rtol=1e-07) | ||
|
|
||
|
|
||
| def test_loss_townsend(): |
There was a problem hiding this comment.
Question for @cwhanse: referring to #1393 (comment), should we have a policy of always testing both array and Series for array-like parameters?
There was a problem hiding this comment.
should we have a policy of always testing both array and Series for array-like parameters?
I think so
There was a problem hiding this comment.
I'd say it's good practice but importance and extent depends on the model and implementation. I'm usually more concerned that we test for compatibility with scalars in "float, array, series" situations because I think it's easier for regressions to slip in with scalars (e.g. by introducing masking).
There was a problem hiding this comment.
This is looking great @abhisheksparikh! I made some suggestions for readability but I'm ok if the group prefers the terser code that more closely matches the reference.
| Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA | ||
|
|
||
| ''' # noqa: E501 | ||
| return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0)) |
There was a problem hiding this comment.
| return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0)) | |
| return np.where(N > 0, 0.5 * S * (1 + 1/N), 0) |
There was a problem hiding this comment.
I'm surprised stickler doesn't complain about this.
| return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0)) | ||
|
|
||
|
|
||
| def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity, |
There was a problem hiding this comment.
"_townsend" or "_townsend_powers"? Often a paper will be cited as "X and Y" if only two authors, "X et al" if 3 or more. Not committed to that pattern here, just providing context. @mikofski?
There was a problem hiding this comment.
It's a tough call, but I believe most folks already refer to this model as the "Townsend snow model"
| Calculates monthly snow loss based on a generalized monthly snow loss model | ||
| discussed in [1]_. |
There was a problem hiding this comment.
"Calculates monthly snow loss based on the Townsend monthly snow loss model [1]_." or Townsend-Powers.
| Snowfall in inches received in a month [in] | ||
|
|
||
| N: numeric | ||
| Number of snowfall events with snowfall > 1" [-] |
| Se_weighted = 1/3 * Se_prev + 2/3 * Se | ||
| gamma = (slant_height * Se_weighted * cosd(surface_tilt)) / \ | ||
| (np.clip((lower_edge_drop_height**2 - Se_weighted**2), a_min=0.01, | ||
| a_max=None) / 2 / tand(angle_of_repose)) |
There was a problem hiding this comment.
In this case, it seems the danger is having negative gamma, which upstream could result in a performance gain (from negative GIT in Eq. 3). To me, clipping gamma to be positive (not zero) is reasonable.
If I understand correctly... I mostly agree, but why can't gamma be 0? Shouldn't that be expected when snow fall is 0? I think mathematically gamma will blow up when both the snow fall is 0 and the drop height goes to 0. But I think in practice python will evaluate 0 / 0 and return nan. Assuming we want to avoid that... Drop height should always be greater than 0, even if just by 1 mm, so maybe that's the point to insert our judgement into the code.
| Se = _townsend_Se(snow_total, snow_events) | ||
| Se_prev = _townsend_Se(snow_total_prev, snow_events_prev) | ||
|
|
||
| Se_weighted = 1/3 * Se_prev + 2/3 * Se |
There was a problem hiding this comment.
| Se = _townsend_Se(snow_total, snow_events) | |
| Se_prev = _townsend_Se(snow_total_prev, snow_events_prev) | |
| Se_weighted = 1/3 * Se_prev + 2/3 * Se | |
| effective_snow = _townsend_Se(snow_total, snow_events) | |
| effective_snow_prev = _townsend_Se(snow_total_prev, snow_events_prev) | |
| effective_snow_weighted = 1/3 * effective_snow_prev + 2/3 * effective_snow |
| gamma = (slant_height * Se_weighted * cosd(surface_tilt)) / \ | ||
| (np.clip((lower_edge_drop_height**2 - Se_weighted**2), a_min=0.01, | ||
| a_max=None) / 2 / tand(angle_of_repose)) |
There was a problem hiding this comment.
| gamma = (slant_height * Se_weighted * cosd(surface_tilt)) / \ | |
| (np.clip((lower_edge_drop_height**2 - Se_weighted**2), a_min=0.01, | |
| a_max=None) / 2 / tand(angle_of_repose)) | |
| drop_height_clipped = np.maximum(lower_edge_drop_height, 0.01) | |
| gamma = ( | |
| slant_height | |
| * effective_snow_weighted | |
| * cosd(surface_tilt)) | |
| / (drop_height_clipped**2 - effective_snow_weighted **2) | |
| * 2 # eqn 5 would benefit from another set of parentheses but I believe | |
| * tand(angle_of_repose) # the 2tan is effectively in the numerator of gamma | |
| ) |
| (np.clip((lower_edge_drop_height**2 - Se_weighted**2), a_min=0.01, | ||
| a_max=None) / 2 / tand(angle_of_repose)) | ||
|
|
||
| GIT = 1 - C2 * np.exp(-gamma) |
There was a problem hiding this comment.
I'd also change this to ground_interference_term and use more lines in the equation below.
* first experimental commit * Added numpy array support * changed var names * changed townsend_Se to private * removed snow.loss_townsend in api.rst * added loss_townsend description in effects_on_pv_system_output.rst * fixed stickler checks * removed rounding of loss and changed to 0-1 range Several other small changes - variable names change, comment change - in response to Kevin's review notes * implementing changes suggested in PR #1251 * removing new line * removing new line * remove new line * Se to effective_snow * adding PR number to whatsnew * address stickler line too long * remove links and noqa E501, and fix long lines * converting to metric system * poa_global from kWh/m2 to Wh/m2 * changing returned loss from kWh to Wh * fixing capacity loss calculation units to keep correct C1 value * convert relative humidity from percent to fraction * neatening docstrings and adjusting variable names * changing eqn 3 percentage loss to fractional loss and adding comment explanation Co-authored-by: Abhishek Parikh <abhishek.parikh2@dnv.com> Co-authored-by: abhisheksparikh <abhi.parikh305@gmail.com>
|
Superseded by #1468 |
docs/sphinx/source/api.rstfor API changes.docs/sphinx/source/whatsnewfor all changes. Includes link to the GitHub Issue with:issue:`num`or this Pull Request with:pull:`num`. Includes contributor name and/or GitHub username (link with:ghuser:`user`).