From 42e0af3a218e2825793c4dae5b385530cb081d9e Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Thu, 4 Sep 2025 13:31:41 -0700 Subject: [PATCH] Switching all calculation functions to be "non-protected" (remove leading underscore) Renamed the omega() and delta() internal calculation functions to be more descriptive and not conflict with the variable names. They are now called "sha" (sunset hour angle) and "declination". --- pyproject.toml | 2 +- src/refet/calcs.py | 130 +++++++++++++++------------------ src/refet/daily.py | 61 +++++++++------- src/refet/hourly.py | 47 ++++++------ tests/test_calcs.py | 174 ++++++++++++++++++++++---------------------- 5 files changed, 206 insertions(+), 208 deletions(-) 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)