diff --git a/pyproject.toml b/pyproject.toml index 67ec3e4..52f6223 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "refet" -version = "0.4.2" +version = "0.5.0" authors = [ { name="Charles Morton", email="charles.morton@dri.edu" }, ] diff --git a/src/refet/calcs.py b/src/refet/calcs.py index 92ec048..3937338 100644 --- a/src/refet/calcs.py +++ b/src/refet/calcs.py @@ -3,7 +3,7 @@ import numpy as np -def _air_pressure(elev, method='asce'): +def air_pressure(elev, method='asce'): """Mean atmospheric pressure at station elevation (Eqs. 3 & 34) Parameters @@ -46,7 +46,7 @@ def _air_pressure(elev, method='asce'): return pair -def _sat_vapor_pressure(temperature): +def sat_vapor_pressure(temperature): """Saturation vapor pressure from temperature (Eq. 7) Parameters @@ -75,7 +75,7 @@ def _sat_vapor_pressure(temperature): return e -def _es_slope(tmean, method='asce'): +def es_slope(tmean, method='asce'): """Slope of the saturation vapor pressure-temperature curve (Eq. 5) Parameters @@ -109,7 +109,7 @@ def _es_slope(tmean, method='asce'): return 4098.0 * 0.6108 * es_slope -def _actual_vapor_pressure(q, pair): +def actual_vapor_pressure(q, pair): """"Actual vapor pressure from specific humidity Parameters @@ -139,7 +139,7 @@ def _actual_vapor_pressure(q, pair): return ea -def _specific_humidity(ea, pair): +def specific_humidity(ea, pair): """"Specific humidity from actual vapor pressure Parameters @@ -169,7 +169,7 @@ def _specific_humidity(ea, pair): return q -def _vpd(es, ea): +def vpd(es, ea): """Vapor pressure deficit Parameters @@ -188,7 +188,7 @@ def _vpd(es, ea): return np.maximum(es - ea, 0) -def _precipitable_water(pair, ea): +def precipitable_water(pair, ea): """Precipitable water in the atmosphere (Eq. D.3) Parameters @@ -207,7 +207,7 @@ def _precipitable_water(pair, ea): return pair * 0.14 * ea + 2.1 -def _doy_fraction(doy): +def doy_fraction(doy): """Fraction of the DOY in the year (Eq. 50) Parameters @@ -224,7 +224,7 @@ def _doy_fraction(doy): return doy * (2 * math.pi / 365) -def _delta(doy, method='asce'): +def declination(doy, method='asce'): """Earth declination (Eq. 51) Parameters @@ -250,12 +250,12 @@ def _delta(doy, method='asce'): """ if method == 'asce': - return 0.409 * np.sin(_doy_fraction(doy) - 1.39) + return 0.409 * np.sin(doy_fraction(doy) - 1.39) elif method == 'refet': return 23.45 * (math.pi / 180) * np.sin(2 * math.pi * (doy + 284) / 365) -def _dr(doy): +def dr(doy): """Inverse square of the Earth-Sun Distance (Eq. 50) Parameters @@ -274,10 +274,10 @@ def _dr(doy): pi * L * d^2 / (ESUN * cos(theta)) -> pi * L / (ESUN * cos(theta) * d) """ - return 1.0 + 0.033 * np.cos(_doy_fraction(doy)) + return 1.0 + 0.033 * np.cos(doy_fraction(doy)) -def _seasonal_correction(doy): +def seasonal_correction(doy): """Seasonal correction for solar time (Eqs. 57 & 58) Parameters @@ -295,7 +295,7 @@ def _seasonal_correction(doy): return 0.1645 * np.sin(2 * b) - 0.1255 * np.cos(b) - 0.0250 * np.sin(b) -def _solar_time_rad(lon, time_mid, sc): +def solar_time_rad(lon, time_mid, sc): """Solar time (i.e. noon is 0) (Eq. 55) Parameters @@ -314,16 +314,16 @@ def _solar_time_rad(lon, time_mid, sc): Notes ----- - This function could be integrated into the _omega() function since they are - always called together (i.e. _omega(_solar_time_rad()). It was built - independently from _omega to eventually support having a separate + This function could be integrated into the sha() function since they are + always called together (i.e. sha(solar_time_rad()). It was built + independently from sha to eventually support having a separate solar_time functions for longitude in degrees. """ return time_mid + (lon * 24 / (2 * math.pi)) + sc - 12 -def _omega(solar_time): +def sha(solar_time): """Solar hour angle (Eq. 55) Parameters @@ -341,10 +341,10 @@ 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) - return _wrap(omega, -math.pi, math.pi) + return wrap(omega, -math.pi, math.pi) -def _wrap(x, x_min, x_max): +def wrap(x, x_min, x_max): """Wrap floating point values into range Parameters @@ -364,7 +364,7 @@ def _wrap(x, x_min, x_max): return np.mod((x - x_min), (x_max - x_min)) + x_min -def _omega_sunset(lat, delta): +def sha_sunset(lat, delta): """Sunset hour angle (Eq. 59) Parameters @@ -383,7 +383,7 @@ def _omega_sunset(lat, delta): return np.arccos(np.clip(-np.tan(lat) * np.tan(delta), -1, 1)) -def _ra_daily(lat, doy, method='asce'): +def ra_daily(lat, doy, method='asce'): """Daily extraterrestrial radiation (Eq. 21) Parameters @@ -408,20 +408,17 @@ def _ra_daily(lat, doy, method='asce'): Equation in Duffie & Beckman (?) uses a solar constant of 1367 W m-2 """ - 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) - ) + delta = declination(doy, method) + omegas = sha_sunset(lat, delta) + theta = omegas * np.sin(lat) * np.sin(delta) + np.cos(lat) * np.cos(delta) * np.sin(omegas) if method == 'asce': - return (24. / math.pi) * 4.92 * _dr(doy) * theta + return (24. / math.pi) * 4.92 * dr(doy) * theta elif method == 'refet': - return (24. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta + return (24. / math.pi) * (1367 * 0.0036) * dr(doy) * theta -def _ra_hourly(lat, lon, doy, time_mid, method='asce'): +def ra_hourly(lat, lon, doy, time_mid, method='asce'): """Hourly extraterrestrial radiation (Eq. 48) Parameters @@ -450,9 +447,9 @@ def _ra_hourly(lat, lon, doy, time_mid, method='asce'): Equation in Duffie & Beckman (?) uses a solar constant of 1367 W m-2 """ - omega = _omega(_solar_time_rad(lon, time_mid, _seasonal_correction(doy))) - delta = _delta(doy, method) - omegas = _omega_sunset(lat, delta) + delta = declination(doy, method) + omega = sha(solar_time_rad(lon, time_mid, seasonal_correction(doy))) + omegas = sha_sunset(lat, delta) # Solar time as start and end of period (Eqs. 53 & 54) # Modify omega1 and omega2 at sunrise and sunset (Eq. 56) @@ -467,12 +464,12 @@ def _ra_hourly(lat, lon, doy, time_mid, method='asce'): ) if method == 'asce': - return (12. / math.pi) * 4.92 * _dr(doy) * theta + return (12. / math.pi) * 4.92 * dr(doy) * theta elif method == 'refet': - return (12. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta + return (12. / math.pi) * (1367 * 0.0036) * dr(doy) * theta -def _rso_daily(ra, ea, pair, doy, lat): +def rso_daily(ra, ea, pair, doy, lat): """Full daily clear sky solar radiation formulation (Appendix D) Parameters @@ -496,19 +493,17 @@ def _rso_daily(ra, ea, pair, doy, lat): """ # 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) + 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) + w = precipitable_water(pair, ea) # Limit sin_beta >= 0.01 so that KB does not go undefined 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 = 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) @@ -517,7 +512,7 @@ def _rso_daily(ra, ea, pair, doy, lat): return (kb + kd) * ra -def _rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): +def rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): """Full hourly clear sky solar radiation formulation (Appendix D) Parameters @@ -540,7 +535,7 @@ def _rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): Calculation method: * 'asce' -- Calculations will follow ASCE-EWRI 2005 [1] equations. * 'refet' -- Calculations will follow RefET software. - Passed through to declination calculation (_delta()). + Passed through to declination calculation (delta()). Returns ------- @@ -548,24 +543,22 @@ def _rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): Hourly clear sky solar radiation [MJ m-2 h-1]. """ - sc = _seasonal_correction(doy) - omega = _omega(_solar_time_rad(lon, time_mid, sc)) + sc = seasonal_correction(doy) + omega = sha(solar_time_rad(lon, time_mid, sc)) # sin of the angle of the sun above the horizon (D.6 and Eq. 62) - delta = _delta(doy, method) + delta = declination(doy, method) sin_beta = np.sin(lat) * np.sin(delta) + np.cos(lat) * np.cos(delta) * np.cos(omega) # Precipitable water - w = _precipitable_water(pair, ea) + w = precipitable_water(pair, ea) # 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 = np.exp( - (-0.00146 * pair) / (kt * sin_beta) - 0.075 * np.power((w / sin_beta), 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) @@ -574,7 +567,7 @@ def _rso_hourly(ra, ea, pair, doy, time_mid, lat, lon, method='asce'): return (kb + kd) * ra -def _rso_simple(ra, elev): +def rso_simple(ra, elev): """Simplified daily/hourly clear sky solar formulation (Eqs. 19 & 45) Parameters @@ -593,7 +586,7 @@ def _rso_simple(ra, elev): return (0.75 + 2E-5 * elev) * ra -def _fcd_daily(rs, rso): +def fcd_daily(rs, rso): """Daytime cloudiness fraction (Eq. 18) Parameters @@ -626,7 +619,7 @@ def _fcd_daily(rs, rso): # return fcd -def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'): +def fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'): """Cloudiness fraction (Eq. 45) Parameters @@ -647,7 +640,7 @@ def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'): Calculation method: * 'asce' -- Calculations will follow ASCE-EWRI 2005 [1] equations. * 'refet' -- Calculations will follow RefET software. - Passed through to declination calculation (_delta()). + Passed through to declination calculation (delta()). Returns ------- @@ -662,12 +655,10 @@ def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'): rso = np.array(rso, copy=True, ndmin=1).astype(np.float64) # DEADBEEF - These values are only needed for identifying low sun angles - sc = _seasonal_correction(doy) - 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) - ) + sc = seasonal_correction(doy) + delta = declination(doy, method) + omega = sha(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)) # As of NumPy 1.7+, ufuncs can take a "where" parameter fcd = np.divide(rs, rso, out=np.ones_like(rs), where=rso != 0) @@ -684,7 +675,7 @@ def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'): return fcd -def _rnl_daily(tmax, tmin, ea, fcd): +def rnl_daily(tmax, tmin, ea, fcd): """Daily net long-wave radiation (Eq. 17) Parameters @@ -710,7 +701,7 @@ def _rnl_daily(tmax, tmin, ea, fcd): ) -def _rnl_hourly(tmean, ea, fcd): +def rnl_hourly(tmean, ea, fcd): """Hourly net long-wave radiation (Eq. 44) Parameters @@ -728,13 +719,10 @@ def _rnl_hourly(tmean, ea, fcd): Hourly net long-wave radiation [MJ m-2 h-1]. """ - return ( - 2.042E-10 * fcd * (0.34 - 0.14 * np.sqrt(ea)) * - np.power((tmean + 273.16), 4) - ) + return 2.042E-10 * fcd * (0.34 - 0.14 * np.sqrt(ea)) * np.power((tmean + 273.16), 4) -def _rn_daily(rs, rnl): +def rn_daily(rs, rnl): """Daily net radiation (Eqs. 15 & 16) Parameters @@ -753,7 +741,7 @@ def _rn_daily(rs, rnl): return 0.77 * rs - rnl -def _rn_hourly(rs, rnl): +def rn_hourly(rs, rnl): """Daily net radiation (Eqs. 42 & 43) Parameters @@ -772,7 +760,7 @@ def _rn_hourly(rs, rnl): return 0.77 * rs - rnl -def _wind_height_adjust(uz, zw): +def wind_height_adjust(uz, zw): """Wind speed at 2 m height based on full logarithmic profile (Eq. 33) Parameters @@ -791,7 +779,7 @@ def _wind_height_adjust(uz, zw): return uz * 4.87 / np.log(67.8 * zw - 5.42) -def _etsz(rn, g, tmean, u2, vpd, es_slope, psy, cn, cd): +def etsz(rn, g, tmean, u2, vpd, es_slope, psy, cn, cd): """Standardized Reference ET [mm] (Eq. 1) Parameters diff --git a/src/refet/daily.py b/src/refet/daily.py index f5af746..0704801 100644 --- a/src/refet/daily.py +++ b/src/refet/daily.py @@ -7,9 +7,22 @@ class Daily(): - def __init__(self, tmin, tmax, rs, uz, zw, elev, lat, doy, ea=None, tdew=None, - 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 @@ -97,13 +110,11 @@ def __init__(self, tmin, tmax, rs, uz, zw, elev, lat, doy, ea=None, tdew=None, # Unit conversions for v, unit in input_units.items(): - setattr( - self, v, units.convert(getattr(self, v), v, unit, timestep='daily') - ) + setattr(self, v, units.convert(getattr(self, v), v, unit, timestep='daily')) # 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) + self.ea = calcs.sat_vapor_pressure(self.tdew) # Rso if rso_type is None: @@ -125,60 +136,56 @@ def __init__(self, tmin, tmax, rs, uz, zw, elev, lat, doy, ea=None, tdew=None, 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) + self.pair = calcs.air_pressure(self.elev, method) # Psychrometric constant (Eq. 4) self.psy = 0.000665 * self.pair # Slope of the saturation vapor pressure-temperature curve - self.es_slope = calcs._es_slope(self.tmean, method) + self.es_slope = calcs.es_slope(self.tmean, method) # Saturation vapor pressure self.es = 0.5 * ( - calcs._sat_vapor_pressure(self.tmax) + - calcs._sat_vapor_pressure(self.tmin) + calcs.sat_vapor_pressure(self.tmax) + + calcs.sat_vapor_pressure(self.tmin) ) # Vapor pressure deficit - self.vpd = calcs._vpd(self.es, self.ea) + self.vpd = calcs.vpd(self.es, self.ea) # Extraterrestrial radiation - self.ra = calcs._ra_daily(self.lat, self.doy, method) + self.ra = calcs.ra_daily(self.lat, self.doy, method) # Clear sky solar radiation # If rso_type is not set, use the method # If rso_type is set, use rso_type directly if rso_type is None : if method.lower() == 'asce': - self.rso = calcs._rso_simple(self.ra, self.elev) + 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.rso = calcs.rso_daily(self.ra, self.ea, self.pair, self.doy, self.lat) elif rso_type.lower() == 'simple': - self.rso = calcs._rso_simple(self.ra, elev) + 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.rso = calcs.rso_daily(self.ra, self.ea, self.pair, self.doy, self.lat) elif rso_type.lower() == 'array': # Use rso array passed to function self.rso = rso # Cloudiness fraction - self.fcd = calcs._fcd_daily(self.rs, self.rso) + self.fcd = calcs.fcd_daily(self.rs, self.rso) # Net long-wave radiation - self.rnl = calcs._rnl_daily(self.tmax, self.tmin, self.ea, self.fcd) + self.rnl = calcs.rnl_daily(self.tmax, self.tmin, self.ea, self.fcd) # Net radiation - self.rn = calcs._rn_daily(self.rs, self.rnl) + self.rn = calcs.rn_daily(self.rs, self.rnl) # Soil heat flux self.g = 0 # Wind speed - self.u2 = calcs._wind_height_adjust(self.uz, self.zw) + self.u2 = calcs.wind_height_adjust(self.uz, self.zw) def etsz(self, surface): """Standardized reference ET @@ -204,7 +211,7 @@ def eto(self): """Grass reference surface""" self.cn = 900 self.cd = 0.34 - return calcs._etsz( + 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 ) @@ -213,7 +220,7 @@ def etr(self): """Alfalfa reference surface""" self.cn = 1600 self.cd = 0.38 - return calcs._etsz( + 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 ) diff --git a/src/refet/hourly.py b/src/refet/hourly.py index dcc97a9..e9891ec 100644 --- a/src/refet/hourly.py +++ b/src/refet/hourly.py @@ -7,9 +7,15 @@ class Hourly(): - def __init__(self, tmean, rs, uz, zw, elev, lat, lon, doy, time, - ea=None, tdew=None, 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]_ @@ -92,13 +98,11 @@ def __init__(self, tmean, rs, uz, zw, elev, lat, lon, doy, time, # Unit conversions for v, unit in input_units.items(): - setattr( - self, v, units.convert(getattr(self, v), v, unit, timestep='hourly') - ) + setattr(self, v, units.convert(getattr(self, v), v, unit, timestep='hourly')) # 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) + 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 @@ -109,31 +113,30 @@ def __init__(self, tmean, rs, uz, zw, elev, lat, lon, doy, time, self.lon *= (math.pi / 180.0) # To match standardized form, psy is calculated from elevation based pair - self.pair = calcs._air_pressure(self.elev, method) + self.pair = calcs.air_pressure(self.elev, method) # Psychrometric constant (Eq. 35) self.psy = 0.000665 * self.pair # Slope of the saturation vapor pressure-temperature curve - self.es_slope = calcs._es_slope(self.tmean, method) + self.es_slope = calcs.es_slope(self.tmean, method) # Saturation vapor pressure - self.es = calcs._sat_vapor_pressure(self.tmean) + self.es = calcs.sat_vapor_pressure(self.tmean) # Vapor pressure deficit self.vpd = self.es - self.ea - # self.vpd = calcs._vpd(self.es, 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': - self.rso = calcs._rso_simple(self.ra, self.elev) + self.rso = calcs.rso_simple(self.ra, self.elev) 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.rso = calcs.rso_hourly( + self.ra, self.ea, self.pair, self.doy, self.time_mid, self.lat, self.lon, method ) # Cloudiness fraction @@ -141,21 +144,21 @@ def __init__(self, tmean, rs, uz, zw, elev, lat, lon, doy, time, # In IN2, "Beta" is computed for the start of the time period, # but "SinBeta" is computed for the midpoint. # Beta (not SinBeta) is used for clamping fcd. - self.fcd = calcs._fcd_hourly( + self.fcd = calcs.fcd_hourly( 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) + self.rnl = calcs.rnl_hourly(self.tmean, self.ea, self.fcd) # Net radiation - self.rn = calcs._rn_hourly(self.rs, self.rnl) + self.rn = calcs.rn_hourly(self.rs, self.rnl) # Soil heat flux (Eqs. 65 and 66) # self.g = self.rn * g_rn # Wind speed - self.u2 = calcs._wind_height_adjust(self.uz, self.zw) + self.u2 = calcs.wind_height_adjust(self.uz, self.zw) def etsz(self, surface): """Standardized reference ET @@ -191,7 +194,7 @@ def eto(self): # Soil heat flux (Eqs. 65 and 66) self.g = self.rn * self.g_rn - return calcs._etsz( + 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 ) @@ -210,7 +213,7 @@ def etr(self): # Soil heat flux (Eqs. 65 and 66) self.g = self.rn * self.g_rn - return calcs._etsz( + 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 ) diff --git a/tests/test_calcs.py b/tests/test_calcs.py index a721fe1..67e85e3 100644 --- a/tests/test_calcs.py +++ b/tests/test_calcs.py @@ -96,13 +96,13 @@ ## Test ancillary functions with positional inputs def test_air_pressure_default(elev=s_args['elev'], pair=s_args['pair_asce']): - assert float(calcs._air_pressure(elev)) == pytest.approx(pair) + 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( @@ -114,7 +114,7 @@ def test_air_pressure_refet(elev=s_args['elev'], pair=s_args['pair']): ] ) def test_sat_vapor_pressure(tdew, ea): - assert float(calcs._sat_vapor_pressure(tdew)) == pytest.approx(ea) + assert float(calcs.sat_vapor_pressure(tdew)) == pytest.approx(ea) @pytest.mark.parametrize( @@ -127,81 +127,77 @@ def test_sat_vapor_pressure(tdew, ea): ] ) def test_specific_humidity_daily(ea, pair, q): - assert float(calcs._specific_humidity(ea, pair)) == pytest.approx(q) + assert float(calcs.specific_humidity(ea, pair)) == pytest.approx(q) -def test_actual_vapor_pressure(q=d_args['q'], pair=s_args['pair'], - ea=d_args['ea']): - assert float(calcs._actual_vapor_pressure(q, pair)) == pytest.approx(ea) +def test_actual_vapor_pressure(q=d_args['q'], pair=s_args['pair'], ea=d_args['ea']): + assert float(calcs.actual_vapor_pressure(q, pair)) == pytest.approx(ea) def test_vpd(es=d_args['es'], ea=d_args['ea']): - assert float(calcs._vpd(es, ea)) == pytest.approx(float(es-ea)) + assert float(calcs.vpd(es, ea)) == pytest.approx(float(es-ea)) # Check that negative VPD's are set to 0 - assert float(calcs._vpd(es, es+1)) == pytest.approx(0) + assert float(calcs.vpd(es, es+1)) == pytest.approx(0) -def test_es_slope_default(tmin=d_args['tmin'], tmax=d_args['tmax'], - es_slope=d_args['es_slope_asce']): - output = calcs._es_slope(0.5 * (tmin + tmax)) +def test_es_slope_default(tmin=d_args['tmin'], tmax=d_args['tmax'], es_slope=d_args['es_slope_asce']): + 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']): - output = calcs._es_slope(0.5 * (tmin + tmax), method='asce') +def test_es_slope_asce(tmin=d_args['tmin'], tmax=d_args['tmax'], es_slope=d_args['es_slope_asce']): + 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']): - output = calcs._es_slope(0.5 * (tmin + tmax), method='refet') +def test_es_slope_refet(tmin=d_args['tmin'], tmax=d_args['tmax'], es_slope=d_args['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) + assert float(calcs.es_slope(tmean, method='asce')) == pytest.approx(es_slope) 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) + 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) + assert float(calcs.precipitable_water(pair, ea)) == pytest.approx(w) def test_doy_fraction(doy=d_args['doy'], expected=d_args['doy_frac']): - assert float(calcs._doy_fraction(doy)) == pytest.approx(expected) + assert float(calcs.doy_fraction(doy)) == pytest.approx(expected) -def test_delta_default(doy=d_args['doy'], delta=d_args['delta_asce']): - assert float(calcs._delta(doy)) == pytest.approx(delta) +def test_declination_default(doy=d_args['doy'], delta=d_args['delta_asce']): + assert float(calcs.declination(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_declination_asce(doy=d_args['doy'], delta=d_args['delta_asce']): + assert float(calcs.declination(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) +def test_declination_refet(doy=d_args['doy'], delta=d_args['delta']): + assert float(calcs.declination(doy, method='refet')) == pytest.approx(delta) def test_dr(doy=d_args['doy'], dr=d_args['dr']): - assert float(calcs._dr(doy)) == pytest.approx(dr) + assert float(calcs.dr(doy)) == pytest.approx(dr) def test_seasonal_correction(doy=d_args['doy'], sc=d_args['sc']): - assert float(calcs._seasonal_correction(doy)) == pytest.approx(sc) + assert float(calcs.seasonal_correction(doy)) == pytest.approx(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']): - assert float(calcs._omega(solar_time)) == pytest.approx(omega) +def test_sha(solar_time=h_args['solar_time'], omega=h_args['omega']): + assert float(calcs.sha(solar_time)) == pytest.approx(omega) @pytest.mark.parametrize( @@ -214,67 +210,64 @@ def test_omega(solar_time=h_args['solar_time'], omega=h_args['omega']): ] ) def test_wrap(x, x_min, x_max, expected): - value = calcs._wrap(x, x_min, x_max) + value = calcs.wrap(x, x_min, x_max) assert pytest.approx(float(value)) == expected -def test_omega_sunset(lat=s_args['lat'], delta=d_args['delta'], - omega_s=d_args['omega_s']): - assert float(calcs._omega_sunset(lat, delta)) == pytest.approx(omega_s) +def test_sha_sunset(lat=s_args['lat'], delta=d_args['delta'], omega_s=d_args['omega_s']): + assert float(calcs.sha_sunset(lat, delta)) == pytest.approx(omega_s) -def test_omega_sunset_high_lat(): - # Test that omega sunset is limited to [0, pi] +def test_sha_sunset_high_lat(): + # Test that sha sunset is limited to [0, pi] # This occurs for angles > ~65 degrees in the summer and ~75 in the winter 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) + assert float(calcs.sha_sunset(lat, calcs.declination(182))) == pytest.approx(math.pi) + assert float(calcs.sha_sunset(-lat, calcs.declination(182))) == pytest.approx(0) + assert float(calcs.sha_sunset(lat, calcs.declination(1))) == pytest.approx(0) 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) + 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) + 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 - assert float(calcs._ra_daily(68 * (math.pi / 180), 1)) == 0 - assert float(calcs._ra_daily(80 * (math.pi / 180), 55)) == 0 + 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'], ra=h_args['ra_asce']): - assert float(calcs._ra_hourly(lat, lon, doy, time)) == pytest.approx(ra) + 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']): - output = calcs._ra_hourly(lat, lon, doy, time, method='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']): - output = calcs._ra_hourly(lat, lon, doy, time, method='refet') + 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 - output = calcs._rso_daily( - calcs._ra_daily(80 * math.pi / 180, 1), ea, pair, 1, 80 * math.pi / 180 - ) + output = calcs.rso_daily(calcs.ra_daily(80 * math.pi / 180, 1), ea, pair, 1, 80 * math.pi / 180) assert float(output) == 0 @@ -282,7 +275,7 @@ 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']): - output = calcs._rso_hourly(ra, ea, pair, doy, time, lat, lon) + output = calcs.rso_hourly(ra, ea, pair, doy, time, lat, lon) assert float(output) == pytest.approx(rso) @@ -290,7 +283,7 @@ 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']): - output = calcs._rso_hourly(ra, ea, pair, doy, time, lat, lon,method='asce') + output = calcs.rso_hourly(ra, ea, pair, doy, time, lat, lon,method='asce') assert float(output) == pytest.approx(rso) @@ -298,7 +291,7 @@ 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']): - output = calcs._rso_hourly(ra, ea, pair, doy, time, lat, lon, method='refet') + output = calcs.rso_hourly(ra, ea, pair, doy, time, lat, lon, method='refet') assert float(output) == pytest.approx(rso) @@ -310,24 +303,24 @@ def test_rso_hourly_refet(ra=h_args['ra'], ea=h_args['ea'], pair=s_args['pair'], ] ) def test_rso_simple(ra, elev, rso): - assert float(calcs._rso_simple(ra, elev)) == pytest.approx(rso) + assert float(calcs.rso_simple(ra, elev)) == pytest.approx(rso) 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) + 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) + 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'], fcd=h_args['fcd_asce']): - output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon) + output = calcs.fcd_hourly(rs, rso, doy, time, lat, lon) assert float(output) == pytest.approx(fcd) @@ -335,7 +328,7 @@ 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']): - output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon, method='asce') + output = calcs.fcd_hourly(rs, rso, doy, time, lat, lon, method='asce') assert float(output) == pytest.approx(fcd) @@ -343,7 +336,7 @@ 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']): - output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon, method='refet') + output = calcs.fcd_hourly(rs, rso, doy, time, lat, lon, method='refet') assert float(output) == pytest.approx(fcd) @@ -351,7 +344,7 @@ 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 - output = calcs._fcd_hourly(rs, rso, doy, time, lat, lon, method='refet') + output = calcs.fcd_hourly(rs, rso, doy, time, lat, lon, method='refet') assert float(output) == pytest.approx(fcd) @@ -362,20 +355,19 @@ def test_fcd_hourly_night(rs=h_args['rs'], rso=h_args['rso'], def test_rnl_daily(tmin=d_args['tmin'], tmax=d_args['tmax'], ea=d_args['ea'], fcd=d_args['fcd'], rnl=d_args['rnl']): - assert float(calcs._rnl_daily(tmax, tmin, ea, fcd)) == pytest.approx(rnl) + assert float(calcs.rnl_daily(tmax, tmin, ea, fcd)) == pytest.approx(rnl) -def test_rnl_hourly(tmean=h_args['tmean'], ea=h_args['ea'], - fcd=h_args['fcd'], rnl=h_args['rnl']): - assert float(calcs._rnl_hourly(tmean, ea, fcd)) == pytest.approx(rnl) +def test_rnl_hourly(tmean=h_args['tmean'], ea=h_args['ea'], fcd=h_args['fcd'], rnl=h_args['rnl']): + assert float(calcs.rnl_hourly(tmean, ea, fcd)) == pytest.approx(rnl) def test_rn_daily(rs=d_args['rs'], rnl=d_args['rnl'], rn=d_args['rn']): - assert float(calcs._rn_daily(rs, rnl)) == pytest.approx(rn) + assert float(calcs.rn_daily(rs, rnl)) == pytest.approx(rn) def test_rn_hourly(rs=h_args['rs'], rnl=h_args['rnl'], rn=h_args['rn']): - assert float(calcs._rn_hourly(rs, rnl)) == pytest.approx(rn) + assert float(calcs.rn_hourly(rs, rnl)) == pytest.approx(rn) @pytest.mark.parametrize( @@ -386,30 +378,38 @@ def test_rn_hourly(rs=h_args['rs'], rnl=h_args['rnl'], rn=h_args['rn']): ] ) def test_wind_height_adjust(uz, zw, u2): - assert float(calcs._wind_height_adjust(uz, zw)) == pytest.approx(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)