diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index 8a8e3ae..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,6 +0,0 @@ -exclude *.txt -include LICENSE.txt README.rst setup.py -recursive-include refet * -recursive-include tests *.py -recursive-exclude tests/data * -exclude MANIFEST.in diff --git a/README.rst b/README.rst index adef99b..31bb7b5 100644 --- a/README.rst +++ b/README.rst @@ -2,7 +2,7 @@ ASCE Standardized Reference Evapotranspiration (ET) =================================================== -|version| |build| |codecov| +|version| |build| NumPy functions for computing daily and hourly reference ET following the ASCE Standardized Reference Evapotranspiration Equations (ASCE2005_). @@ -21,16 +21,18 @@ The raw input data is available `here C - ea = 0.6108 * math.exp(17.27 * tdew_c / (tdew_c + 237.3)) # kPa - # ea = refet.calcs._saturated_vapor_pressure(tdew_c) + # The actual vapor pressure could be computed from the dew point temperature below + # or the tdew can be passed directly to the function + # Convert the dew point temperature to Celsius + # tdew = units._f2c(49.84) + # ea = 0.6108 * math.exp(17.27 * tdew / (tdew + 237.3)) + # ea = refet.calcs._sat_vapor_pressure(tdew) etr = refet.Daily( - tmin=66.65, tmax=102.80, ea=ea, rs=674.07, uz=4.80, + tmin=66.65, tmax=102.80, tdew=49.84, rs=674.07, uz=4.80, zw=3, elev=1208.5, lat=39.4575, doy=182, method='asce', - input_units={'tmin': 'F', 'tmax': 'F', 'rs': 'Langleys', 'uz': 'mph', - 'lat': 'deg'} + input_units={'tmin': 'F', 'tmax': 'F', 'tdew': 'F', 'rs': 'Langleys', + 'uz': 'mph', 'lat': 'deg'} ).etr() print(f'ETr: {float(etr):.2f} mm') @@ -64,7 +66,6 @@ Required Parameters (hourly & daily) ======== ========== ==================================================== Variable Type Description [default units] ======== ========== ==================================================== -ea ndarray Actual vapor pressure [kPa] uz ndarray Wind speed [m s-1] zw float Wind speed height [m] elev ndarray Elevation [m] @@ -72,6 +73,18 @@ lat ndarray Latitude [degrees] doy ndarray Day of year ======== ========== ==================================================== +Required Ea Parameters (hourly & daily) +--------------------------------------------------- + +Either the "ea" or "tdew" parameter must be set + +======== ========== ==================================================== +Variable Type Description [default units] +======== ========== ==================================================== +ea ndarray Actual vapor pressure [kPa] +tdew ndarray Dew point temperature [C] +======== ========== ==================================================== + Required Daily Parameters ------------------------- @@ -185,6 +198,3 @@ References .. |version| image:: https://badge.fury.io/py/refet.svg :alt: Latest version on PyPI :target: https://badge.fury.io/py/refet -.. |codecov| image:: https://codecov.io/gh/WSWUP/refet/branch/main/graphs/badge.svg - :alt: Coverage Status - :target: https://codecov.io/gh/WSWUP/refet diff --git a/pyproject.toml b/pyproject.toml index 4e87998..ded321d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "refet" -version = "0.3.12" +version = "0.4.0" authors = [ { name="Charles Morton", email="charles.morton@dri.edu" }, ] diff --git a/src/refet/calcs.py b/src/refet/calcs.py index c1ce55d..92ec048 100644 --- a/src/refet/calcs.py +++ b/src/refet/calcs.py @@ -17,7 +17,7 @@ def _air_pressure(elev, method='asce'): Returns ------- - pair : ndarray + ndarray Air pressure [kPa]. Notes @@ -42,6 +42,7 @@ def _air_pressure(elev, method='asce'): np.power(pair, 9.8 / (0.0065 * 286.9), out=pair) # np.power(pair, 5.26, out=pair) pair *= 101.3 + return pair @@ -55,7 +56,7 @@ def _sat_vapor_pressure(temperature): Returns ------- - e : ndarray + ndarray Saturation vapor pressure [kPa]. Notes @@ -70,6 +71,7 @@ def _sat_vapor_pressure(temperature): e *= 17.27 np.exp(e, out=e) e *= 0.6108 + return e @@ -87,22 +89,24 @@ def _es_slope(tmean, method='asce'): Returns ------- - es_slope : ndarray + ndarray Slope [kPa C-1]. Notes ----- - 4098 * 0.6108 * exp(17.27 * T / (T + 237.3)) / ((T + 237.3) ** 2)) + The difference between the two methods is in the rounding of the coefficient. + + Intentionally (partially) recomputing es from tmean instead of taking es as + function parameter since this es is computed directly from tmean instead of + as the average of the min and max temperature saturation vapor pressures. """ - if method == 'refet': - es_slope = ( - 4098.0 * _sat_vapor_pressure(tmean) / np.power(tmean + 237.3, 2)) - elif method == 'asce': - es_slope = ( - 2503.0 * np.exp(17.27 * tmean / (tmean + 237.3)) / - np.power(tmean + 237.3, 2)) - return es_slope + es_slope = np.exp(17.27 * tmean / (tmean + 237.3)) / np.power(tmean + 237.3, 2) + + if method == 'asce': + return 2503.0 * es_slope + elif method == 'refet': + return 4098.0 * 0.6108 * es_slope def _actual_vapor_pressure(q, pair): @@ -117,7 +121,7 @@ def _actual_vapor_pressure(q, pair): Returns ------- - ea : ndarray + ndarray Actual vapor pressure [kPa]. Notes @@ -131,6 +135,7 @@ def _actual_vapor_pressure(q, pair): np.reciprocal(ea, out=ea) ea *= pair ea *= q + return ea @@ -146,7 +151,7 @@ def _specific_humidity(ea, pair): Returns ------- - q : ndarray + ndarray Specific humidity [kg/kg]. Notes @@ -160,6 +165,7 @@ def _specific_humidity(ea, pair): np.reciprocal(q, out=q) q *= ea q *= 0.622 + return q @@ -179,7 +185,6 @@ def _vpd(es, ea): Vapor pressure deficit [kPa]. """ - return np.maximum(es - ea, 0) @@ -216,7 +221,7 @@ def _doy_fraction(doy): DOY fraction [radians]. """ - return doy * (2.0 * math.pi / 365) + return doy * (2 * math.pi / 365) def _delta(doy, method='asce'): @@ -286,7 +291,7 @@ def _seasonal_correction(doy): Seasonal correction [hour] """ - b = 2 * math.pi * (doy - 81.) / 364. + b = 2 * math.pi * (doy - 81) / 364.0 return 0.1645 * np.sin(2 * b) - 0.1255 * np.cos(b) - 0.0250 * np.sin(b) @@ -328,7 +333,7 @@ def _omega(solar_time): Returns ------- - omega : ndarray + ndarray Hour angle [radians]. """ @@ -336,8 +341,7 @@ def _omega(solar_time): # Need to adjust omega so that the values go from -pi to pi # Values outside this range are wrapped (i.e. -3*pi/2 -> pi/2) - omega = _wrap(omega, -math.pi, math.pi) - return omega + return _wrap(omega, -math.pi, math.pi) def _wrap(x, x_min, x_max): @@ -395,7 +399,7 @@ def _ra_daily(lat, doy, method='asce'): Returns ------- - ra : ndarray + ndarray Daily extraterrestrial radiation [MJ m-2 d-1]. Notes @@ -406,14 +410,15 @@ def _ra_daily(lat, doy, method='asce'): """ delta = _delta(doy, method) omegas = _omega_sunset(lat, delta) - theta = (omegas * np.sin(lat) * np.sin(delta) + - np.cos(lat) * np.cos(delta) * np.sin(omegas)) + theta = ( + omegas * np.sin(lat) * np.sin(delta) + + np.cos(lat) * np.cos(delta) * np.sin(omegas) + ) if method == 'asce': - ra = (24. / math.pi) * 4.92 * _dr(doy) * theta + return (24. / math.pi) * 4.92 * _dr(doy) * theta elif method == 'refet': - ra = (24. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta - return ra + return (24. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta def _ra_hourly(lat, lon, doy, time_mid, method='asce'): @@ -436,7 +441,7 @@ def _ra_hourly(lat, lon, doy, time_mid, method='asce'): Returns ------- - ra : ndarray + ndarray Hourly extraterrestrial radiation [MJ m-2 h-1]. Notes @@ -458,12 +463,13 @@ def _ra_hourly(lat, lon, doy, time_mid, method='asce'): # Extraterrestrial radiation (Eq. 48) theta = ( ((omega2 - omega1) * np.sin(lat) * np.sin(delta)) + - (np.cos(lat) * np.cos(delta) * (np.sin(omega2) - np.sin(omega1)))) + (np.cos(lat) * np.cos(delta) * (np.sin(omega2) - np.sin(omega1))) + ) + if method == 'asce': - ra = (12. / math.pi) * 4.92 * _dr(doy) * theta + return (12. / math.pi) * 4.92 * _dr(doy) * theta elif method == 'refet': - ra = (12. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta - return ra + return (12. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta def _rso_daily(ra, ea, pair, doy, lat): @@ -484,29 +490,31 @@ def _rso_daily(ra, ea, pair, doy, lat): Returns ------- - rso : ndarray + ndarray Daily clear sky solar radiation [MJ m-2 d-1] """ # sin of the angle of the sun above the horizon (D.5 and Eq. 62) sin_beta_24 = np.sin( - 0.85 + 0.3 * lat * np.sin(_doy_fraction(doy) - 1.39) - - 0.42 * np.power(lat, 2)) - sin_beta_24 = np.maximum(sin_beta_24, 0.1) + 0.85 + 0.3 * lat * np.sin(_doy_fraction(doy) - 1.39) - 0.42 * np.power(lat, 2) + ) # Precipitable water w = _precipitable_water(pair, ea) - # Clearness index for direct beam radiation (Eq. D.2) # Limit sin_beta >= 0.01 so that KB does not go undefined - kb = (0.98 * np.exp((-0.00146 * pair) / sin_beta_24 - - 0.075 * np.power((w / sin_beta_24), 0.4))) + sin_beta_24 = np.maximum(sin_beta_24, 0.1) + + # Clearness index for direct beam radiation (Eq. D.2) + kb = np.exp( + (-0.00146 * pair) / sin_beta_24 - 0.075 * np.power((w / sin_beta_24), 0.4) + ) + kb *= 0.98 # Transmissivity index for diffuse radiation (Eq. D.4) kd = np.minimum(-0.36 * kb + 0.35, 0.82 * kb + 0.18) - rso = ra * (kb + kd) - return rso + return (kb + kd) * ra def _rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): @@ -536,7 +544,7 @@ def _rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): Returns ------- - rso : ndarray + ndarray Hourly clear sky solar radiation [MJ m-2 h-1]. """ @@ -545,25 +553,25 @@ def _rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): # sin of the angle of the sun above the horizon (D.6 and Eq. 62) delta = _delta(doy, method) - sin_beta = ( - np.sin(lat) * np.sin(delta) + - np.cos(lat) * np.cos(delta) * np.cos(omega)) + sin_beta = np.sin(lat) * np.sin(delta) + np.cos(lat) * np.cos(delta) * np.cos(omega) # Precipitable water w = _precipitable_water(pair, ea) - # Clearness index for direct beam radiation (Eq. D.2) # Limit sin_beta >= 0.01 so that KB does not go undefined + sin_beta = np.maximum(sin_beta, 0.01) + + # Clearness index for direct beam radiation (Eq. D.2) kt = 1.0 - kb = 0.98 * np.exp( - (-0.00146 * pair) / (kt * np.maximum(sin_beta, 0.01)) - - 0.075 * np.power((w / np.maximum(sin_beta, 0.01)), 0.4)) + kb = np.exp( + (-0.00146 * pair) / (kt * sin_beta) - 0.075 * np.power((w / sin_beta), 0.4) + ) + kb *= 0.98 # Transmissivity index for diffuse radiation (Eq. D.4) kd = np.minimum(-0.36 * kb + 0.35, 0.82 * kb + 0.18) - rso = ra * (kb + kd) - return rso + return (kb + kd) * ra def _rso_simple(ra, elev): @@ -578,12 +586,11 @@ def _rso_simple(ra, elev): Returns ------- - rso : ndarray + ndarray Clear sky solar radiation [MJ m-2 d-1 or MJ m-2 h-1]. """ - rso = (0.75 + 2E-5 * elev) * ra - return rso + return (0.75 + 2E-5 * elev) * ra def _fcd_daily(rs, rso): @@ -659,8 +666,8 @@ def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'): delta = _delta(doy, method) omega = _omega(_solar_time_rad(lon, time_mid, sc)) beta = np.arcsin( - np.sin(lat) * np.sin(delta) + - np.cos(lat) * np.cos(delta) * np.cos(omega)) + np.sin(lat) * np.sin(delta) + np.cos(lat) * np.cos(delta) * np.cos(omega) + ) # As of NumPy 1.7+, ufuncs can take a "where" parameter fcd = np.divide(rs, rso, out=np.ones_like(rs), where=rso != 0) @@ -699,7 +706,8 @@ def _rnl_daily(tmax, tmin, ea, fcd): """ return ( 4.901E-9 * fcd * (0.34 - 0.14 * np.sqrt(ea)) * - 0.5 * (np.power(tmax + 273.16, 4) + np.power(tmin + 273.16, 4))) + 0.5 * (np.power(tmax + 273.16, 4) + np.power(tmin + 273.16, 4)) + ) def _rnl_hourly(tmean, ea, fcd): @@ -722,7 +730,8 @@ def _rnl_hourly(tmean, ea, fcd): """ return ( 2.042E-10 * fcd * (0.34 - 0.14 * np.sqrt(ea)) * - np.power((tmean + 273.16), 4)) + np.power((tmean + 273.16), 4) + ) def _rn_daily(rs, rnl): @@ -813,4 +822,5 @@ def _etsz(rn, g, tmean, u2, vpd, es_slope, psy, cn, cd): """ return ( (0.408 * es_slope * (rn - g) + (psy * cn * u2 * vpd / (tmean + 273))) / - (es_slope + psy * (cd * u2 + 1))) + (es_slope + psy * (cd * u2 + 1)) + ) diff --git a/src/refet/daily.py b/src/refet/daily.py index 27df419..f5af746 100644 --- a/src/refet/daily.py +++ b/src/refet/daily.py @@ -7,8 +7,9 @@ class Daily(): - def __init__(self, tmin, tmax, ea, rs, uz, zw, elev, lat, doy, - method='asce', rso_type=None, rso=None, input_units={}): + def __init__(self, tmin, tmax, rs, uz, zw, elev, lat, doy, ea=None, tdew=None, + method='asce', rso_type=None, rso=None, input_units={}, + ): """ASCE Daily Standardized Reference Evapotranspiration (ET) Arguments @@ -17,8 +18,6 @@ def __init__(self, tmin, tmax, ea, rs, uz, zw, elev, lat, doy, Minimum daily temperature [C]. tmax : ndarray Maximum daily temperature [C]. - ea : ndarray - Actual vapor pressure [kPa]. rs : ndarray Incoming shortwave solar radiation [MJ m-2 day-1]. uz : ndarray @@ -31,6 +30,10 @@ def __init__(self, tmin, tmax, ea, rs, uz, zw, elev, lat, doy, Latitude [degrees]. doy : ndarray Day of year. + ea : ndarray, optional + Actual vapor pressure [kPa]. Either ea or tdew must be set. + tdew : ndarray, optional + Mean daily dew point temperature [C]. method : {'asce' (default), 'refet'}, optional Specifies which calculation method to use. * 'asce' -- Calculations will follow ASCE-EWRI 2005 [1]_ equations. @@ -68,11 +71,12 @@ def __init__(self, tmin, tmax, ea, rs, uz, zw, elev, lat, doy, Evapotranspiration Task Committee Rep., ASCE Reston, Va. """ + if method.lower() not in ['asce', 'refet']: + raise ValueError('method must be "asce" or "refet"') # Convert all inputs to NumPy arrays self.tmin = np.array(tmin, copy=True, ndmin=1) self.tmax = np.array(tmax, copy=True, ndmin=1) - self.ea = np.array(ea, copy=True, ndmin=1) self.rs = np.array(rs, copy=True, ndmin=1) self.uz = np.array(uz, copy=True, ndmin=1) self.elev = np.array(elev, copy=True, ndmin=1) @@ -80,16 +84,28 @@ def __init__(self, tmin, tmax, ea, rs, uz, zw, elev, lat, doy, self.zw = zw self.doy = doy - # Unit conversion + # Use Ea directly if it is set, otherwise try to compute from Tdew + if ea is not None: + self.ea = np.array(ea, copy=True, ndmin=1) + self.tdew = None + elif tdew is not None: + self.tdew = np.array(tdew, copy=True, ndmin=1) + self.ea = None + else: + # TODO: Check if there is a better exception to raise + raise Exception('Either "ea" or "tdew" parameter must be set') + + # Unit conversions for v, unit in input_units.items(): setattr( - self, v, - units.convert(getattr(self, v), v, unit, timestep='daily') + self, v, units.convert(getattr(self, v), v, unit, timestep='daily') ) - if method.lower() not in ['asce', 'refet']: - raise ValueError('method must be "asce" or "refet"') + # Compute Ea after handling unit conversions so that Tdew is in Celsius + if self.ea is None and self.tdew is not None: + self.ea = calcs._sat_vapor_pressure(self.tdew) + # Rso if rso_type is None: pass elif rso_type.lower() not in ['simple', 'full', 'array']: @@ -105,19 +121,23 @@ def __init__(self, tmin, tmax, ea, rs, uz, zw, elev, lat, doy, # expecting the default units to be radians self.lat *= (math.pi / 180.0) + # Mean daily air temperature + self.tmean = 0.5 * (self.tmax + self.tmin) + # To match standardized form, pair is calculated from elevation self.pair = calcs._air_pressure(self.elev, method) # Psychrometric constant (Eq. 4) self.psy = 0.000665 * self.pair - self.tmean = 0.5 * (self.tmax + self.tmin) + # Slope of the saturation vapor pressure-temperature curve self.es_slope = calcs._es_slope(self.tmean, method) - # Saturated vapor pressure + # Saturation vapor pressure self.es = 0.5 * ( calcs._sat_vapor_pressure(self.tmax) + - calcs._sat_vapor_pressure(self.tmin)) + calcs._sat_vapor_pressure(self.tmin) + ) # Vapor pressure deficit self.vpd = calcs._vpd(self.es, self.ea) @@ -133,12 +153,14 @@ def __init__(self, tmin, tmax, ea, rs, uz, zw, elev, lat, doy, self.rso = calcs._rso_simple(self.ra, self.elev) elif method.lower() == 'refet': self.rso = calcs._rso_daily( - self.ra, self.ea, self.pair, self.doy, self.lat) + self.ra, self.ea, self.pair, self.doy, self.lat + ) elif rso_type.lower() == 'simple': self.rso = calcs._rso_simple(self.ra, elev) elif rso_type.lower() == 'full': self.rso = calcs._rso_daily( - self.ra, self.ea, self.pair, self.doy, self.lat) + self.ra, self.ea, self.pair, self.doy, self.lat + ) elif rso_type.lower() == 'array': # Use rso array passed to function self.rso = rso @@ -184,7 +206,8 @@ def eto(self): self.cd = 0.34 return calcs._etsz( rn=self.rn, g=self.g, tmean=self.tmean, u2=self.u2, vpd=self.vpd, - es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd) + es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd + ) def etr(self): """Alfalfa reference surface""" @@ -192,4 +215,5 @@ def etr(self): self.cd = 0.38 return calcs._etsz( rn=self.rn, g=self.g, tmean=self.tmean, u2=self.u2, vpd=self.vpd, - es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd) + es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd + ) diff --git a/src/refet/hourly.py b/src/refet/hourly.py index 246a51d..2ec3c57 100644 --- a/src/refet/hourly.py +++ b/src/refet/hourly.py @@ -7,8 +7,9 @@ class Hourly(): - def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, - method='asce', input_units={}): + def __init__(self, tmean, rs, uz, zw, elev, lat, lon, doy, time, + ea=None, tdew=None, method='asce', input_units={}, + ): """ASCE Hourly Standardized Reference Evapotranspiration (ET) .. warning:: Cloudiness fraction at night is not being computed per [1]_ @@ -17,8 +18,6 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, --------- tmean : ndarray Average hourly temperature [C]. - ea : ndarray - Actual vapor pressure [kPa]. rs : ndarray Shortwave solar radiation [MJ m-2 hr-1]. uz : ndarray @@ -35,6 +34,10 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, Day of year. time : ndarray UTC hour at start of time period. + ea : ndarray, optional + Actual vapor pressure [kPa]. Either ea or tdew parameter must be set. + tdew : ndarray, optional + Average hourly dew point temperature [C]. method : {'asce' (default), 'refet'}, optional Specifies which calculation method to use. * 'asce' -- Calculations will follow ASCE-EWRI 2005 [1]_ equations. @@ -60,10 +63,11 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, Evapotranspiration Task Committee Rep., ASCE Reston, Va. """ + if method.lower() not in ['asce', 'refet']: + raise ValueError('method must be "asce" or "refet"') # Convert all inputs to NumPy arrays self.tmean = np.array(tmean, copy=True, ndmin=1) - self.ea = np.array(ea, copy=True, ndmin=1) self.rs = np.array(rs, copy=True, ndmin=1) self.uz = np.array(uz, copy=True, ndmin=1) self.elev = np.array(elev, copy=True, ndmin=1) @@ -75,15 +79,27 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, self.zw = zw self.doy = doy - # Unit conversion + # Use Ea directly if it is set, otherwise try to compute from Tdew + if ea is not None: + self.ea = np.array(ea, copy=True, ndmin=1) + self.tdew = None + elif tdew is not None: + self.tdew = np.array(tdew, copy=True, ndmin=1) + self.ea = None + else: + # TODO: Check if there is a better exception to raise + raise Exception('Either "ea" or "tdew" parameter must be set') + + # Unit conversions for v, unit in input_units.items(): + print(v, unit) setattr( - self, v, - units.convert(getattr(self, v), v, unit, timestep='hourly') + self, v, units.convert(getattr(self, v), v, unit, timestep='hourly') ) - if method.lower() not in ['asce', 'refet']: - raise ValueError('method must be "asce" or "refet"') + # Compute Ea after handling unit conversions so that Tdew is in Celsius + if self.ea is None and self.tdew is not None: + self.ea = calcs._sat_vapor_pressure(self.tdew) # The input angles are converted to degrees by default in units.convert # They need to be converted back to radians for the calc functions @@ -99,16 +115,18 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, # Psychrometric constant (Eq. 35) self.psy = 0.000665 * self.pair - self.es = calcs._sat_vapor_pressure(self.tmean) + # Slope of the saturation vapor pressure-temperature curve self.es_slope = calcs._es_slope(self.tmean, method) + # Saturation vapor pressure + self.es = calcs._sat_vapor_pressure(self.tmean) + # Vapor pressure deficit self.vpd = self.es - self.ea # self.vpd = calcs._vpd(self.es, ea) # Extraterrestrial radiation - self.ra = calcs._ra_hourly( - self.lat, self.lon, self.doy, self.time_mid, method) + self.ra = calcs._ra_hourly(self.lat, self.lon, self.doy, self.time_mid, method) # Clear sky solar radiation if method == 'asce': @@ -116,7 +134,8 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, elif method == 'refet': self.rso = calcs._rso_hourly( self.ra, self.ea, self.pair, self.doy, self.time_mid, self.lat, - self.lon, method) + self.lon, method + ) # Cloudiness fraction # Intentionally not using time_mid to match Beta value in IN2 file @@ -124,7 +143,8 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, # but "SinBeta" is computed for the midpoint. # Beta (not SinBeta) is used for clamping fcd. self.fcd = calcs._fcd_hourly( - self.rs, self.rso, self.doy, self.time, self.lat, self.lon, method) + self.rs, self.rso, self.doy, self.time, self.lat, self.lon, method + ) # Net long-wave radiation self.rnl = calcs._rnl_hourly(self.tmean, self.ea, self.fcd) @@ -174,7 +194,8 @@ def eto(self): return calcs._etsz( rn=self.rn, g=self.g, tmean=self.tmean, u2=self.u2, vpd=self.vpd, - es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd) + es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd + ) def etr(self): """Tall (alfalfa) reference surface""" @@ -192,4 +213,5 @@ def etr(self): return calcs._etsz( rn=self.rn, g=self.g, tmean=self.tmean, u2=self.u2, vpd=self.vpd, - es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd) + es_slope=self.es_slope, psy=self.psy, cn=self.cn, cd=self.cd + ) diff --git a/src/refet/units.py b/src/refet/units.py index 5f12ef2..0e6a6e8 100644 --- a/src/refet/units.py +++ b/src/refet/units.py @@ -52,7 +52,7 @@ def convert(values, variable, unit, timestep=None): # Convert input values to expected units # TODO: Split these into separate functions - if variable in ['tmean', 'tmin', 'tmax']: + if variable in ['tmean', 'tmin', 'tmax', 'tdew']: if unit.lower() in ['f', 'fahrenheit']: values -= 32 values *= (5.0 / 9) diff --git a/tests/test_calcs.py b/tests/test_calcs.py index 97e47d3..a721fe1 100644 --- a/tests/test_calcs.py +++ b/tests/test_calcs.py @@ -99,29 +99,33 @@ def test_air_pressure_default(elev=s_args['elev'], pair=s_args['pair_asce']): assert float(calcs._air_pressure(elev)) == pytest.approx(pair) def test_air_pressure_asce(elev=s_args['elev'], pair=s_args['pair_asce']): - assert float(calcs._air_pressure( - elev, method='asce')) == pytest.approx(pair) + assert float(calcs._air_pressure(elev, method='asce')) == pytest.approx(pair) def test_air_pressure_refet(elev=s_args['elev'], pair=s_args['pair']): - assert float(calcs._air_pressure( - elev, method='refet')) == pytest.approx(pair) + assert float(calcs._air_pressure(elev, method='refet')) == pytest.approx(pair) @pytest.mark.parametrize( 'tdew, ea', - [[d_args['tdew'], d_args['ea']], - [h_args['tmean'], h_args['es']], - [h_args['tdew'], h_args['ea']]]) + [ + [d_args['tdew'], d_args['ea']], + [h_args['tmean'], h_args['es']], + [h_args['tdew'], h_args['ea']], + ] +) def test_sat_vapor_pressure(tdew, ea): assert float(calcs._sat_vapor_pressure(tdew)) == pytest.approx(ea) @pytest.mark.parametrize( 'ea, pair, q', - [[d_args['ea'], s_args['pair'], d_args['q']], - [d_args['ea'], s_args['pair_asce'], d_args['q_asce']], - [h_args['ea'], s_args['pair'], h_args['q']], - [h_args['ea'], s_args['pair_asce'], h_args['q_asce']]]) + [ + [d_args['ea'], s_args['pair'], d_args['q']], + [d_args['ea'], s_args['pair_asce'], d_args['q_asce']], + [h_args['ea'], s_args['pair'], h_args['q']], + [h_args['ea'], s_args['pair_asce'], h_args['q_asce']], + ] +) def test_specific_humidity_daily(ea, pair, q): assert float(calcs._specific_humidity(ea, pair)) == pytest.approx(q) @@ -139,22 +143,31 @@ def test_vpd(es=d_args['es'], ea=d_args['ea']): def test_es_slope_default(tmin=d_args['tmin'], tmax=d_args['tmax'], es_slope=d_args['es_slope_asce']): - assert float(calcs._es_slope( - 0.5 * (tmin + tmax))) == pytest.approx(es_slope) + output = calcs._es_slope(0.5 * (tmin + tmax)) + assert float(output) == pytest.approx(es_slope) + def test_es_slope_asce(tmin=d_args['tmin'], tmax=d_args['tmax'], es_slope=d_args['es_slope_asce']): - assert float(calcs._es_slope( - 0.5 * (tmin + tmax), method='asce')) == pytest.approx(es_slope) + output = calcs._es_slope(0.5 * (tmin + tmax), method='asce') + assert float(output) == pytest.approx(es_slope) + def test_es_slope_refet(tmin=d_args['tmin'], tmax=d_args['tmax'], es_slope=d_args['es_slope']): - assert float(calcs._es_slope( - 0.5 * (tmin + tmax), method='refet')) == pytest.approx(es_slope) + output = calcs._es_slope(0.5 * (tmin + tmax), method='refet') + assert float(output) == pytest.approx(es_slope) + + +def test_es_slope_hourly_asce(tmean=h_args['tmean'], es_slope=h_args['es_slope_asce']): + assert float(calcs._es_slope(tmean, method='asce')) == pytest.approx(es_slope) -def test_precipitable_water(pair=s_args['pair'], ea=d_args['ea'], - w=d_args['w']): +def test_es_slope_hourly_refet(tmean=h_args['tmean'], es_slope=h_args['es_slope']): + assert float(calcs._es_slope(tmean, method='refet')) == pytest.approx(es_slope) + + +def test_precipitable_water(pair=s_args['pair'], ea=d_args['ea'], w=d_args['w']): assert float(calcs._precipitable_water(pair, ea)) == pytest.approx(w) @@ -165,9 +178,11 @@ def test_doy_fraction(doy=d_args['doy'], expected=d_args['doy_frac']): def test_delta_default(doy=d_args['doy'], delta=d_args['delta_asce']): assert float(calcs._delta(doy)) == pytest.approx(delta) + def test_delta_asce(doy=d_args['doy'], delta=d_args['delta_asce']): assert float(calcs._delta(doy, method='asce')) == pytest.approx(delta) + def test_delta_refet(doy=d_args['doy'], delta=d_args['delta']): assert float(calcs._delta(doy, method='refet')) == pytest.approx(delta) @@ -182,8 +197,7 @@ def test_seasonal_correction(doy=d_args['doy'], sc=d_args['sc']): def test_solar_time_rad(lon=s_args['lon'], time_mid=h_args['time'], sc=d_args['sc'], expected=h_args['solar_time']): - assert float(calcs._solar_time_rad( - lon, time_mid, sc)) == pytest.approx(expected) + assert float(calcs._solar_time_rad(lon, time_mid, sc)) == pytest.approx(expected) def test_omega(solar_time=h_args['solar_time'], omega=h_args['omega']): @@ -192,10 +206,12 @@ def test_omega(solar_time=h_args['solar_time'], omega=h_args['omega']): @pytest.mark.parametrize( 'x, x_min, x_max, expected', - [[1.2, 1.2, 1.5, 1.2], - [1.1, 1.2, 1.5, 1.4], - [1.6, 1.2, 1.5, 1.3], - [2.0, 1.2, 1.5, 1.4]] + [ + [1.2, 1.2, 1.5, 1.2], + [1.1, 1.2, 1.5, 1.4], + [1.6, 1.2, 1.5, 1.3], + [2.0, 1.2, 1.5, 1.4], + ] ) def test_wrap(x, x_min, x_max, expected): value = calcs._wrap(x, x_min, x_max) @@ -216,18 +232,17 @@ def test_omega_sunset_high_lat(): 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'], - ra=d_args['ra_asce']): +def test_ra_daily_default(lat=s_args['lat'], doy=d_args['doy'], ra=d_args['ra_asce']): assert float(calcs._ra_daily(lat, doy)) == pytest.approx(ra) -def test_ra_daily_asce(lat=s_args['lat'], doy=d_args['doy'], - ra=d_args['ra_asce']): - assert float(calcs._ra_daily( - lat, doy, method='asce')) == pytest.approx(ra) + +def test_ra_daily_asce(lat=s_args['lat'], doy=d_args['doy'], ra=d_args['ra_asce']): + assert float(calcs._ra_daily(lat, doy, method='asce')) == pytest.approx(ra) + 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) + 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 @@ -235,66 +250,65 @@ def test_ra_daily_zero(): 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'], - ra=h_args['ra_asce']): - assert float(calcs._ra_hourly( - lat, lon, doy, time)) == pytest.approx(ra) +def test_ra_hourly_default(lat=s_args['lat'], lon=s_args['lon'], doy=h_args['doy'], + time=h_args['time_mid'], ra=h_args['ra_asce']): + assert float(calcs._ra_hourly(lat, lon, doy, time)) == pytest.approx(ra) -def test_ra_hourly_asce(lat=s_args['lat'], lon=s_args['lon'], - doy=h_args['doy'], time=h_args['time_mid'], - ra=h_args['ra_asce']): - assert float(calcs._ra_hourly( - lat, lon, doy, time, method='asce')) == pytest.approx(ra) +def test_ra_hourly_asce(lat=s_args['lat'], lon=s_args['lon'], doy=h_args['doy'], + time=h_args['time_mid'], ra=h_args['ra_asce']): + output = calcs._ra_hourly(lat, lon, doy, time, method='asce') + assert float(output) == pytest.approx(ra) -def test_ra_hourly_refet(lat=s_args['lat'], lon=s_args['lon'], - doy=h_args['doy'], time=h_args['time_mid'], - ra=h_args['ra']): - assert float(calcs._ra_hourly( - lat, lon, doy, time, method='refet')) == pytest.approx(ra) +def test_ra_hourly_refet(lat=s_args['lat'], lon=s_args['lon'], doy=h_args['doy'], + time=h_args['time_mid'], ra=h_args['ra']): + output = calcs._ra_hourly(lat, lon, doy, time, method='refet') + assert float(output) == pytest.approx(ra) def test_rso_daily(ra=d_args['ra'], ea=d_args['ea'], pair=s_args['pair'], doy=d_args['doy'], lat=s_args['lat'], rso=d_args['rso']): - assert float(calcs._rso_daily( - ra, ea, pair, doy, lat)) == pytest.approx(rso) + assert float(calcs._rso_daily(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 + output = calcs._rso_daily( + calcs._ra_daily(80 * math.pi / 180, 1), ea, pair, 1, 80 * math.pi / 180 + ) + assert float(output) == 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'], rso=h_args['rso_asce']): - assert float(calcs._rso_hourly( - ra, ea, pair, doy, time, lat, lon)) == pytest.approx(rso) + output = calcs._rso_hourly(ra, ea, pair, doy, time, lat, lon) + assert float(output) == pytest.approx(rso) + def test_rso_hourly_asce(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'], rso=h_args['rso_asce']): - assert float(calcs._rso_hourly( - ra, ea, pair, doy, time, lat, lon, - method='asce')) == pytest.approx(rso) + output = calcs._rso_hourly(ra, ea, pair, doy, time, lat, lon,method='asce') + assert float(output) == pytest.approx(rso) + def test_rso_hourly_refet(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'], rso=h_args['rso']): - assert float(calcs._rso_hourly( - ra, ea, pair, doy, time, lat, lon, - method='refet')) == pytest.approx(rso) + output = calcs._rso_hourly(ra, ea, pair, doy, time, lat, lon, method='refet') + assert float(output) == pytest.approx(rso) @pytest.mark.parametrize( 'ra, elev, rso', - [[d_args['ra'], s_args['elev'], d_args['rso_simple']], - [h_args['ra'], s_args['elev'], h_args['rso_simple']]]) + [ + [d_args['ra'], s_args['elev'], d_args['rso_simple']], + [h_args['ra'], s_args['elev'], h_args['rso_simple']], + ] +) def test_rso_simple(ra, elev, rso): assert float(calcs._rso_simple(ra, elev)) == pytest.approx(rso) @@ -313,29 +327,32 @@ 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'], fcd=h_args['fcd_asce']): - assert float(calcs._fcd_hourly( - rs, rso, doy, time, lat, lon)) == pytest.approx(fcd) + output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon) + assert float(output) == pytest.approx(fcd) + def test_fcd_hourly_asce(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'], fcd=h_args['fcd_asce']): - assert float(calcs._fcd_hourly( - rs, rso, doy, time, lat, lon, method='asce')) == pytest.approx(fcd) + output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon, method='asce') + assert float(output) == pytest.approx(fcd) + def test_fcd_hourly_refet(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'], fcd=h_args['fcd']): - assert float(calcs._fcd_hourly( - rs, rso, doy, time, lat, lon, method='refet')) == pytest.approx(fcd) + output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon, method='refet') + assert float(output) == pytest.approx(fcd) + def test_fcd_hourly_night(rs=h_args['rs'], rso=h_args['rso'], doy=h_args['doy'], time=6, lat=s_args['lat'], lon=s_args['lon'], fcd=1): # For now, check that nighttime fcd values are set to 1 - assert float(calcs._fcd_hourly( - rs, rso, doy, time, lat, lon, method='refet')) == pytest.approx(fcd) + output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon, method='refet') + assert float(output) == pytest.approx(fcd) # Test function when rs/rso have different shapes than lat/lon/doy/time @@ -363,33 +380,36 @@ def test_rn_hourly(rs=h_args['rs'], rnl=h_args['rnl'], rn=h_args['rn']): @pytest.mark.parametrize( 'uz, zw, u2', - [[d_args['uz'], s_args['zw'], d_args['u2']], - [h_args['uz'], s_args['zw'], h_args['u2']]]) + [ + [d_args['uz'], s_args['zw'], d_args['u2']], + [h_args['uz'], s_args['zw'], h_args['u2']], + ] +) def test_wind_height_adjust(uz, zw, u2): assert float(calcs._wind_height_adjust(uz, zw)) == pytest.approx(u2) def test_wind_height_adjust_2m(uz=2.5, zw=2.0, u2=2.5): - assert float(calcs._wind_height_adjust( - uz, zw)) == pytest.approx(u2, abs=0.001) + assert float(calcs._wind_height_adjust(uz, zw)) == pytest.approx(u2, abs=0.001) @pytest.mark.parametrize( 'rn, g, tmean, u2, es, ea, es_slope, pair, cn, cd, etsz', - [[d_args['rn'], 0, 0.5 * (d_args['tmin'] + d_args['tmax']), d_args['u2'], - d_args['es'], d_args['ea'], d_args['es_slope'], s_args['pair'], - 900, 0.34, d_args['eto']], - [d_args['rn'], 0, 0.5 * (d_args['tmin'] + d_args['tmax']), d_args['u2'], - d_args['es'], d_args['ea'], d_args['es_slope'], s_args['pair'], - 1600, 0.38, d_args['etr']], - [h_args['rn'], 0.1 * h_args['rn'], h_args['tmean'], h_args['u2'], - h_args['es'], h_args['ea'], h_args['es_slope'], s_args['pair'], - 37.0, 0.24, h_args['eto']], - [h_args['rn'], 0.04 * h_args['rn'], h_args['tmean'], h_args['u2'], - h_args['es'], h_args['ea'], h_args['es_slope'], s_args['pair'], - 66.0, 0.25, h_args['etr']], - ]) + [ + [d_args['rn'], 0, 0.5 * (d_args['tmin'] + d_args['tmax']), d_args['u2'], + d_args['es'], d_args['ea'], d_args['es_slope'], s_args['pair'], + 900, 0.34, d_args['eto']], + [d_args['rn'], 0, 0.5 * (d_args['tmin'] + d_args['tmax']), d_args['u2'], + d_args['es'], d_args['ea'], d_args['es_slope'], s_args['pair'], + 1600, 0.38, d_args['etr']], + [h_args['rn'], 0.1 * h_args['rn'], h_args['tmean'], h_args['u2'], + h_args['es'], h_args['ea'], h_args['es_slope'], s_args['pair'], + 37.0, 0.24, h_args['eto']], + [h_args['rn'], 0.04 * h_args['rn'], h_args['tmean'], h_args['u2'], + h_args['es'], h_args['ea'], h_args['es_slope'], s_args['pair'], + 66.0, 0.25, h_args['etr']], + ] +) def test_etsz(rn, g, tmean, u2, es, ea, es_slope, pair, cn, cd, etsz): - output = calcs._etsz(rn, g, tmean, u2, es - ea, es_slope, 0.000665 * pair, - cn, cd) + output = calcs._etsz(rn, g, tmean, u2, es - ea, es_slope, 0.000665 * pair, cn, cd) assert float(output) == pytest.approx(etsz, abs=0.001) diff --git a/tests/test_daily.py b/tests/test_daily.py index 7877083..90ce8bf 100644 --- a/tests/test_daily.py +++ b/tests/test_daily.py @@ -38,9 +38,9 @@ # Test full daily functions with positional inputs def test_refet_daily_input_positions(): etr = Daily( - d_args['tmin'], d_args['tmax'], d_args['ea'], d_args['rs'], + d_args['tmin'], d_args['tmax'], d_args['rs'], d_args['uz'], s_args['zw'], s_args['elev'], s_args['lat'], - d_args['doy'], 'asce', + d_args['doy'], d_args['ea'], None, 'asce', ).etr() assert float(etr) == pytest.approx(d_args['etr_asce']) @@ -241,3 +241,37 @@ def test_refet_daily_lat_rad(): input_units={'lat': 'rad'}, ).etr() assert float(etr) == pytest.approx(d_args['etr_asce']) + + +def test_refet_daily_tdew(): + """Check that Ea can be computed from tdew""" + etr = Daily( + tmin=d_args['tmin'], tmax=d_args['tmax'], tdew=d_args['tdew'], + rs=d_args['rs'], uz=d_args['uz'], zw=s_args['zw'], + elev=s_args['elev'], lat=s_args['lat'], doy=d_args['doy'], + method='asce', + ).etr() + assert float(etr) == pytest.approx(d_args['etr_asce']) + + +def test_refet_daily_ea_tdew_exception(): + """Check that an exception is raised if ea and tdew are not set""" + with pytest.raises(Exception): + Daily( + tmin=d_args['tmin'], tmax=d_args['tmax'], + rs=d_args['rs'], uz=d_args['uz'], zw=s_args['zw'], + elev=s_args['elev'], lat=s_args['lat'], doy=d_args['doy'], + method='asce', + ) + + +def test_refet_daily_tdew(): + """Check that Tdew is converted to C before computing Ea""" + etr = Daily( + tmin=d_args['tmin'], tmax=d_args['tmax'], + tdew=units._c2f(d_args['tdew']), + rs=d_args['rs'], uz=d_args['uz'], zw=s_args['zw'], + elev=s_args['elev'], lat=s_args['lat'], doy=d_args['doy'], + method='asce', input_units={'tdew': 'F'}, + ).etr() + assert float(etr) == pytest.approx(d_args['etr_asce']) \ No newline at end of file diff --git a/tests/test_hourly.py b/tests/test_hourly.py index 98277be..6c0a4d0 100644 --- a/tests/test_hourly.py +++ b/tests/test_hourly.py @@ -39,9 +39,9 @@ # Test full hourly functions with positional inputs def test_refet_hourly_input_positions(): etr = Hourly( - h_args['tmean'], h_args['ea'], h_args['rs'], h_args['uz'], + h_args['tmean'], h_args['rs'], h_args['uz'], s_args['zw'], s_args['elev'], s_args['lat'], s_args['lon'], - h_args['doy'], h_args['time'], 'asce', + h_args['doy'], h_args['time'], h_args['ea'], None, 'asce', ).etr() assert float(etr) == pytest.approx(h_args['etr_asce']) @@ -216,3 +216,37 @@ def test_refet_hourly_lon_rad(): input_units={'lon': 'rad'}, ).etr() assert float(etr) == pytest.approx(h_args['etr_asce']) + + +def test_refet_hourly_tdew(): + """Check that Ea can be computed from tdew""" + etr = Hourly( + tmean=h_args['tmean'], tdew=h_args['tdew'], rs=h_args['rs'], + uz=h_args['uz'], zw=s_args['zw'], elev=s_args['elev'], + lat=s_args['lat'], lon=s_args['lon'], doy=h_args['doy'], + time=h_args['time'], + ).etr() + assert float(etr) == pytest.approx(h_args['etr_asce']) + + +def test_refet_hourly_ea_tdew_exception(): + """Check that an exception is raised if ea and tdew are not set""" + with pytest.raises(Exception): + Hourly( + tmean=h_args['tmean'], rs=h_args['rs'], + uz=h_args['uz'], zw=s_args['zw'], elev=s_args['elev'], + lat=s_args['lat'], lon=s_args['lon'], doy=h_args['doy'], + time=h_args['time'], + ) + + +def test_refet_hourly_tdew(): + """Check that Tdew is converted to C before computing Ea""" + etr = Hourly( + tmean=h_args['tmean'], rs=h_args['rs'], + tdew=units._c2f(h_args['tdew']), + uz=h_args['uz'], zw=s_args['zw'], elev=s_args['elev'], + lat=s_args['lat'], lon=s_args['lon'], doy=h_args['doy'], + time=h_args['time'], input_units={'tdew': 'F'}, + ).etr() + assert float(etr) == pytest.approx(h_args['etr_asce'])