Fixed error when Rso is zero (for high latitudes and/or winter DOY)

Fcd will be set to 1 when Rso is 0.  This mimics the approach that is in the hourly calculation.

refet/__init__.py

@@ -1,4 +1,4 @@
 from .daily import Daily
 from .hourly import Hourly
 
-__version__ = "0.3.8"
+__version__ = "0.3.9"

refet/calcs.py

@@ -600,8 +600,23 @@ def _fcd_daily(rs, rso):
  -------
  ndarray
 
+ Notes
+ -----
+ fcd will be set to 1 when Rso is 0
+
  """
- return 1.35 * np.clip(rs / rso, 0.3, 1.0) - 0.35
+ rs = np.array(rs, copy=True, ndmin=1).astype(np.float64)
+ rso = np.array(rso, copy=True, ndmin=1).astype(np.float64)
+
+ # As of NumPy 1.7+, ufuncs can take a "where" parameter
+ fcd = np.divide(rs, rso, out=np.ones_like(rs), where=rso != 0)
+ return 1.35 * np.clip(fcd, 0.3, 1.0) - 0.35
+ # return 1.35 * np.clip(rs / rso, 0.3, 1.0) - 0.35
+
+ # # DEADBEEF
+ # fcd = np.ones(rso.shape)
+ # fcd[rso > 0] = 1.35 * np.clip(rs[rso > 0] / rso[rso > 0], 0.3, 1) - 0.35
+ # return fcd
 
 
 def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'):
@@ -629,7 +644,11 @@ def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'):
 
  Returns
  -------
- fcd : ndarray
+ ndarray
+
+ Notes
+ -----
+ fcd will be set to 1 when Rso is 0
 
  """
  rs = np.array(rs, copy=True, ndmin=1).astype(np.float64)
@@ -643,10 +662,14 @@ def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'):
  np.sin(lat) * np.sin(delta) +
  np.cos(lat) * np.cos(delta) * np.cos(omega))
 
- fcd = np.ones(rso.shape)
- fcd[rso > 0] = 1.35 * np.clip(rs[rso > 0] / rso[rso > 0], 0.3, 1) - 0.35
+ # As of NumPy 1.7+, ufuncs can take a "where" parameter
+ fcd = np.divide(rs, rso, out=np.ones_like(rs), where=rso != 0)
+ fcd = 1.35 * np.clip(fcd, 0.3, 1.0) - 0.35
+ # DEADBEEF
+ # fcd = np.ones(rso.shape)
+ # fcd[rso > 0] = 1.35 * np.clip(rs[rso > 0] / rso[rso > 0], 0.3, 1) - 0.35
 
- # For now just set fcd to 1 for low sun angles
+ # For now set fcd to 1 for low sun angles also
  # DEADBEEF - Still need to get daytime value of fcd when beta > 0.3
  # Get closest value in time (array space) when beta > 0.3
  fcd[beta < 0.3] = 1

tests/test_calcs.py

@@ -1,5 +1,6 @@
 import math
 
+import numpy as np
 import pytest
 
 import refet.calcs as calcs
@@ -204,11 +205,13 @@ def test_omega_sunset(lat=s_args['lat'], delta=d_args['delta'],
  assert float(calcs._omega_sunset(lat, delta)) == pytest.approx(omega_s)
 
 
-def test_omega_sunset_high_lat(lat=70*(math.pi/180), delta=d_args['delta']):
+def test_omega_sunset_high_lat():
  # Test that omega sunset is limited to [0, pi]
  # This occurs for angles > ~65 degrees in the summer and ~75 in the winter
- assert float(calcs._omega_sunset(lat, delta)) == pytest.approx(math.pi)
- assert float(calcs._omega_sunset(-lat, delta)) == pytest.approx(0)
+ lat = 70 * (math.pi / 180)
+ assert float(calcs._omega_sunset(lat, calcs._delta(182))) == pytest.approx(math.pi)
+ assert float(calcs._omega_sunset(-lat, calcs._delta(182))) == pytest.approx(0)
+ assert float(calcs._omega_sunset(lat, calcs._delta(1))) == pytest.approx(0)
 
 
 def test_ra_daily_default(lat=s_args['lat'], doy=d_args['doy'],
@@ -224,6 +227,11 @@ def test_ra_daily_refet(lat=s_args['lat'], doy=d_args['doy'], ra=d_args['ra']):
  assert float(calcs._ra_daily(
  lat, doy, method='refet')) == pytest.approx(ra)
 
+def test_ra_daily_zero():
+ # Ra can go to zero for winter DOY and/or high latitudes
+ assert float(calcs._ra_daily(68 * (math.pi / 180), 1)) == 0
+ assert float(calcs._ra_daily(80 * (math.pi / 180), 55)) == 0
+
 
 def test_ra_hourly_default(lat=s_args['lat'], lon=s_args['lon'],
  doy=h_args['doy'], time=h_args['time_mid'],
@@ -250,6 +258,13 @@ def test_rso_daily(ra=d_args['ra'], ea=d_args['ea'], pair=s_args['pair'],
  ra, ea, pair, doy, lat)) == pytest.approx(rso)
 
 
+def test_rso_daily_ra_zero(ea=d_args['ea'], pair=s_args['pair']):
+ # Rso can go to zero for winter DOY and/or high latitudes when Ra is zero
+ rso = calcs._rso_daily(calcs._ra_daily(80 * math.pi / 180, 1), ea, pair, 1,
+ 80 * math.pi / 180)
+ assert float(rso) == 0
+
+
 def test_rso_hourly_default(ra=h_args['ra'], ea=h_args['ea'], pair=s_args['pair'],
  doy=h_args['doy'], time=h_args['time_mid'],
  lat=s_args['lat'], lon=s_args['lon'],
@@ -286,6 +301,12 @@ def test_fcd_daily(rs=d_args['rs'], rso=d_args['rso'], fcd=d_args['fcd']):
  assert float(calcs._fcd_daily(rs, rso)) == pytest.approx(fcd)
 
 
+def test_fcd_daily_rso_zero():
+ # Check that fcd returns 1 when Rso is zero
+ assert float(calcs._fcd_daily(1, 0)) == pytest.approx(1.0)
+ assert float(calcs._fcd_daily(0, 0)) == pytest.approx(1.0)
+
+
 def test_fcd_hourly_default(rs=h_args['rs'], rso=h_args['rso'],
  doy=h_args['doy'], time=h_args['time_mid'],
  lat=s_args['lat'], lon=s_args['lon'],