Townsend snow #1251

abhisheksparikh · 2021-07-02T16:19:16Z

Closes Implement monthly snow loss model #1246
I am familiar with the contributing guidelines
Tests added
Updates entries to docs/sphinx/source/api.rst for API changes.
Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for 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`).
New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
Pull request is nearly complete and ready for detailed review.
Maintainer: Appropriate GitHub Labels and Milestone are assigned to the Pull Request and linked Issue.

cwhanse

Thanks @abhisheksparikh this is a great start. Some comments for you to consider as you work through the formatting.

cwhanse · 2021-07-02T17:53:31Z

pvlib/snow.py

+    '''
+    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):


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.

pvlib/snow.py

abhisheksparikh · 2022-03-15T16:03:42Z

Excuse me for the long delay! I will be working on resolving the issues of this PR in the coming couple of days.

into townsend_snow

kandersolar

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?

pvlib/snow.py

Several other small changes - variable names change, comment change - in response to Kevin's review notes

kandersolar

More nitpicking. I think after this it will be close!

kandersolar · 2022-03-18T17:59:29Z

docs/sphinx/source/whatsnew/v0.9.1.rst

@@ -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`


Suggested change

* Added Townsend Powers Snow loss model in :py:func:`pvlib.snow`

* Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.townsend`

kandersolar · 2022-03-18T18:05:45Z

pvlib/snow.py

+    -------
+    loss : numeric
+        Average monthly DC capacity loss due to snow coverage [%]
+


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.

kandersolar · 2022-03-18T18:09:35Z

pvlib/snow.py

+    lower_edge_drop_height : float
+        Drop height from array edge to ground [in]
+
+    P : float


This P slipped through the cracks :)

kandersolar · 2022-03-18T18:11:18Z

pvlib/snow.py

+
+    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


Good catch that temperature is in Kelvin 👍

kandersolar · 2022-03-18T18:12:11Z

pvlib/snow.py

+        Relative humidity [%]
+
+    temp_air : numeric
+        Ambient temperature [°C]


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?

kandersolar · 2022-03-18T18:28:00Z

pvlib/snow.py

+        Row length in the slanted plane of array dimension [in]
+
+    lower_edge_drop_height : float
+        Drop height from array edge to ground [in]


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.

I prefer metric system over everything else :)

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?

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.

I'm OK asking users to do some unit conversions to keep consistent units for similar functions in pvlib

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.

kandersolar · 2022-03-18T18:28:48Z

pvlib/snow.py

+    Returns
+    -------
+    loss : numeric
+        Average monthly DC capacity loss due to snow coverage [%]


Suggested change

Average monthly DC capacity loss due to snow coverage [%]

Average monthly DC capacity loss fraction due to snow coverage

kandersolar · 2022-03-18T18:33:49Z

pvlib/snow.py

+    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))


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.

I think a_min is set at 0.01 to prevent gamma from blowing up.

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 :)

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

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,

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?

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.

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.

kandersolar · 2022-03-18T18:35:44Z

pvlib/snow.py

+
+    Parameters
+    ----------
+    snow_total : numeric


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

kandersolar · 2022-03-18T18:39:09Z

pvlib/tests/test_snow.py

+    np.testing.assert_allclose(expected, actual, rtol=1e-07)
+
+
+def test_loss_townsend():


Question for @cwhanse: referring to #1393 (comment), should we have a policy of always testing both array and Series for array-like parameters?

should we have a policy of always testing both array and Series for array-like parameters?

I think so

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).

wholmgren

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.

wholmgren · 2022-03-22T17:01:53Z

pvlib/snow.py

+       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))


Suggested change

return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0))

return np.where(N > 0, 0.5 * S * (1 + 1/N), 0)

I'm surprised stickler doesn't complain about this.

wholmgren · 2022-03-22T17:07:20Z

pvlib/snow.py

+    return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0))
+
+
+def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,


"_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?

It's a tough call, but I believe most folks already refer to this model as the "Townsend snow model"

wholmgren · 2022-03-22T17:09:46Z

pvlib/snow.py

+    Calculates monthly snow loss based on a generalized monthly snow loss model
+    discussed in [1]_.


"Calculates monthly snow loss based on the Townsend monthly snow loss model [1]_." or Townsend-Powers.

wholmgren · 2022-03-22T17:13:50Z

pvlib/snow.py

+        Snowfall in inches received in a month [in]
+
+    N: numeric
+        Number of snowfall events with snowfall > 1" [-]


1" here too

wholmgren · 2022-03-22T17:35:59Z

pvlib/snow.py

+    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))


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.

wholmgren · 2022-03-22T17:51:30Z

pvlib/snow.py

+    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


Suggested change

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

wholmgren · 2022-03-22T17:58:15Z

pvlib/snow.py

+    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))


Suggested change

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

)

wholmgren · 2022-03-22T17:59:37Z

pvlib/snow.py

+        (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)


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>

kandersolar · 2022-09-12T22:45:47Z

Superseded by #1468

Abhishek Parikh added 2 commits June 24, 2021 21:19

first experimental commit

0b736d4

Added numpy array support

3f2c40d

cwhanse reviewed Jul 2, 2021

View reviewed changes

changed var names

974d516

mikofski added the enhancement label Jul 19, 2021

mikofski added this to the 0.9.0 milestone Jul 19, 2021

mikofski assigned abhisheksparikh Jul 19, 2021

changed townsend_Se to private

655fa79

wholmgren modified the milestones: 0.9.0, 0.9.1 Aug 7, 2021

cwhanse mentioned this pull request Mar 14, 2022

v0.9.1 release planning #1424

Closed

14 tasks

abhisheksparikh added 4 commits March 16, 2022 21:35

removed snow.loss_townsend in api.rst

c0c8b4d

Merge branch 'master' of https://github.com/abhisheksparikh/pvlib-python

f0e1f3a

into townsend_snow

added loss_townsend description in effects_on_pv_system_output.rst

78716f7

fixed stickler checks

5e06a40

kandersolar reviewed Mar 17, 2022

View reviewed changes

abhisheksparikh added 2 commits March 18, 2022 00:34

removed rounding of loss and changed to 0-1 range

a5c5315

Several other small changes - variable names change, comment change - in response to Kevin's review notes

Merge branch 'master' into townsend_snow

f79e906

kandersolar reviewed Mar 18, 2022

View reviewed changes

wholmgren reviewed Mar 22, 2022

View reviewed changes

kandersolar modified the milestones: 0.9.1, 0.9.2 Mar 29, 2022

reepoi added a commit to reepoi/pvlib-python that referenced this pull request May 26, 2022

implementing changes suggested in PR pvlib#1251

c2d6d4f

reepoi mentioned this pull request Jun 7, 2022

Townsend snow (building on #1251) #1468

Merged

8 tasks

kandersolar modified the milestones: 0.9.2, 0.9.3 Aug 19, 2022

kandersolar closed this Sep 12, 2022

kandersolar mentioned this pull request Jan 17, 2023

Corrections to Townsend snow model #1636

Closed

	* Added Townsend Powers Snow loss model in :py:func:`pvlib.snow`
	* Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.townsend`

	Average monthly DC capacity loss due to snow coverage [%]
	Average monthly DC capacity loss fraction due to snow coverage

		np.testing.assert_allclose(expected, actual, rtol=1e-07)


		def test_loss_townsend():

	return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0))
	return np.where(N > 0, 0.5 * S * (1 + 1/N), 0)

		return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0))


		def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,

		Calculates monthly snow loss based on a generalized monthly snow loss model
		discussed in [1]_.

Townsend snow #1251

Townsend snow #1251

Uh oh!

Conversation

abhisheksparikh commented Jul 2, 2021 • edited by mikofski Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

cwhanse left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

abhisheksparikh commented Mar 15, 2022

Uh oh!

kandersolar left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

kandersolar left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

abhisheksparikh Mar 22, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

abhisheksparikh Mar 22, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mikofski Mar 22, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wholmgren Mar 22, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

abhisheksparikh commented Jul 2, 2021 •

edited by mikofski

Loading

abhisheksparikh Mar 22, 2022 •

edited

Loading

abhisheksparikh Mar 22, 2022 •

edited

Loading

mikofski Mar 22, 2022 •

edited

Loading

wholmgren Mar 22, 2022 •

edited

Loading