Hello Im using a rasperry 3 weewx v3.9.2 and i want to add some special data and also use this for some pin outputs (heating on-of for melting snow and so one).
Therefore i installed your filepille extension and add some code. I now have to import data like outTemp to calculate an wetbulb temperatur. I tried many parses and it worked in wxformulas but only for HTML but there are no recordet data for images or max min values. Does anybody know the right import and parsing to put in filepille.py, for a value like current.outTemp? -- You received this message because you are subscribed to the Google Groups "weewx-user" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/weewx-user/a5b7be97-4dc6-4db1-8cb1-38043aed4812%40googlegroups.com.
# # Copyright (c) 2019 Tom Keffer <[email protected]> # # See the file LICENSE.txt for your full rights. # """ A WeeWX that parses a file, adding its contents to a WeeWX record. To use: 1. Include a stanza in your weewx.conf configuration file: [FilePile] filename = /var/tmp/filepile.txt unit_system = METRIC # Or, 'US' or 'METRICWX' # Map from incoming names, to WeeWX names. [[label_map]] temp1 = extraTemp1 humid1 = extraHumid1 2. Add the FilePile service to the list of data_services to be run: [Engine] [[Services]] ... data_services = user.filepile.FilePile 3. Put this file (filepile.py) in your WeeWX user subdirectory. For example, if you installed using setup.py, cp filepile.py /home/weewx/bin/user 4. Have your external data source write values to the file ('/var/tmp/filepile.txt' in the example above) using the following format: key = value where 'key' is an observation name, and 'value' is its value. The value can be 'None'. """ #ad imports (more guessing than knowing) import syslog import weewx import weewx.units import Adafruit_DHT import weewx.wxformulas import weewx.uwxutils from weewx.wxengine import StdService from weeutil.weeutil import to_float from weewx.units import INHG_PER_MBAR, METER_PER_FOOT, METER_PER_MILE, MM_PER_INCH from weewx.units import CtoK, CtoF, FtoC class FilePile(StdService): """WeeWX service for augmenting a record with data parsed from a file.""" def __init__(self, engine, config_dict): # Initialize my superclass: super(FilePile, self).__init__(engine, config_dict) # Extract our stanza from the configuration dicdtionary filepile_dict = config_dict.get('FilePile', {}) # Get the location of the file ... self.filename = filepile_dict.get('filename', '/var/tmp/filepile.txt') # ... and the unit system it will use unit_system_name = filepile_dict.get('unit_system', 'METRICWX').strip().upper() # Make sure we know about the unit system. If not, raise an exception. if unit_system_name not in weewx.units.unit_constants: raise ValueError("FilePile: Unknown unit system: %s" % unit_system_name) # Use the numeric code for the unit system self.unit_system = weewx.units.unit_constants[unit_system_name] # Mapping from variable names to weewx names self.label_map = filepile_dict.get('label_map', {}) syslog.syslog(syslog.LOG_INFO, "filepile: Using %s with the '%s' unit system" % (self.filename, unit_system_name)) syslog.syslog(syslog.LOG_INFO, "filepile: Label map is %s" % self.label_map) # Bind to the NEW_ARCHIVE_RECORD event self.bind(weewx.NEW_ARCHIVE_RECORD, self.new_archive_record) def new_archive_record(self, event): """ #Add New code # Code with formola for wettbulb. Addet to wxformula (works with $current.cloudbase) def wetbulp(t_C, rh): #Calculate the wet bulp temperatur #t_C - temperature in degrees Celsius #rh - relative humidity [0-100] dp_ce = dewpointC(t_C, rh) if dp_ce is None: return None #wbt2 = t_C * math.atan(0.151977 * math.pow((rh + 8.313659), 0.5)) + math.atan(t_C + rh) - math.atan(rh - 1.676331) + 0.00391838 * math.pow(rh, 1.5) * math.atan(0.023101 * rh) - 4.686035 wbt = -5.809 + 0.058 * rh + 0.697 * t_C + 0.003 * rh * t_C return wbt """ sensorcam = Adafruit_DHT.DHT22 pin = 26 humiditycam, temperaturecam = Adafruit_DHT.read_retry(sensorcam, pin) #Not working frases for extraTemp_1: #weewx.wxformulas.wetbulp #weewx.wxformulas.wetbulp(t_C, rh) #weewx.current.outTemp #weewx.observation type.outTemp extraTemp_1 = None extraTemp_2 = round(temperaturecam, 2) extraTemp_3 = None extraHumid_1 = round(humiditycam, 2) snow_ = None hail_ = None if not isinstance(extraTemp_1 , (float, int)): extraTemp_1 = 3 if not isinstance(extraTemp_2 , (float, int)): extraTemp_2 = None elif not isinstance(extraTemp_3 , (float, int)): extraTemp_3 = None elif not isinstance(extraHumid_1 , (float, int)): extraHumid_1 = None elif not isinstance(snow_ , (float, int)): snow_ = None elif not isinstance(hail_ , (float, int)): hail_ = None meineextras = open(self.filename, 'w') meineextras.write('extraTemp1 = ' + str(extraTemp_1) + '\nextraTemp2 = ' + str(extraTemp_2) + '\nextraTemp3 = ' + str(extraTemp_3) + '\nextraHumid1 = ' + str(extraHumid_1) + '\nsnow = ' + str(snow_) + '\nhail = ' + str(hail_)) meineextras.close() '''Original Code''' new_record_data = {} try: with open(self.filename, 'r') as fd: for line in fd: eq_index = line.find('=') # Ignore all lines that do not have an equal sign if eq_index == -1: continue name = line[:eq_index].strip() value = line[eq_index + 1:].strip() new_record_data[self.label_map.get(name, name)] = to_float(value) # Supply a unit system if one wasn't included in the file if 'usUnits' not in new_record_data: new_record_data['usUnits'] = self.unit_system # Convert the new values to the same unit system as the record target_data = weewx.units.to_std_system(new_record_data, event.record['usUnits']) # Add the converted values to the record: event.record.update(target_data) except IOError as e: syslog.syslog(syslog.LOG_ERR, "FilePile: Cannot open file. Reason: %s" % e)
# # Copyright (c) 2009-2015 Tom Keffer <[email protected]> # # See the file LICENSE.txt for your full rights. # """Various weather related formulas and utilities.""" import math import syslog import time import weewx.uwxutils from weewx.units import INHG_PER_MBAR, METER_PER_FOOT, METER_PER_MILE, MM_PER_INCH from weewx.units import CtoK, CtoF, FtoC def dewpointF(T, R): """Calculate dew point. T: Temperature in Fahrenheit R: Relative humidity in percent. Returns: Dewpoint in Fahrenheit Examples: >>> print "%.1f" % dewpointF(68, 50) 48.7 >>> print "%.1f" % dewpointF(32, 50) 15.5 >>> print "%.1f" % dewpointF(-10, 50) -23.5 """ if T is None or R is None: return None TdC = dewpointC(FtoC(T), R) return CtoF(TdC) if TdC is not None else None def dewpointC(T, R): """Calculate dew point. http://en.wikipedia.org/wiki/Dew_point T: Temperature in Celsius R: Relative humidity in percent. Returns: Dewpoint in Celsius """ if T is None or R is None: return None R = R / 100.0 try: _gamma = 17.27 * T / (237.7 + T) + math.log(R) TdC = 237.7 * _gamma / (17.27 - _gamma) except (ValueError, OverflowError): TdC = None return TdC def windchillF(T_F, V_mph): """Calculate wind chill. http://www.nws.noaa.gov/om/cold/wind_chill.shtml T_F: Temperature in Fahrenheit V_mph: Wind speed in mph Returns Wind Chill in Fahrenheit """ if T_F is None or V_mph is None: return None # only valid for temperatures below 50F and wind speeds over 3.0 mph if T_F >= 50.0 or V_mph <= 3.0: return T_F WcF = 35.74 + 0.6215 * T_F + (-35.75 + 0.4275 * T_F) * math.pow(V_mph, 0.16) return WcF def windchillC(T_C, V_kph): """Wind chill, metric version. T: Temperature in Celsius V: Wind speed in kph Returns wind chill in Celsius""" if T_C is None or V_kph is None: return None T_F = CtoF(T_C) V_mph = 0.621371192 * V_kph WcF = windchillF(T_F, V_mph) return FtoC(WcF) if WcF is not None else None def heatindexF(T, R): """Calculate heat index. http://www.nws.noaa.gov/om/heat/heat_index.shtml T: Temperature in Fahrenheit R: Relative humidity in percent Returns heat index in Fahrenheit Examples: >>> print heatindexF(75.0, 50.0) 75.0 >>> print heatindexF(80.0, 50.0) 80.8029049 >>> print heatindexF(80.0, 95.0) 86.3980618 >>> print heatindexF(90.0, 50.0) 94.5969412 >>> print heatindexF(90.0, 95.0) 126.6232036 """ if T is None or R is None: return None # Formula only valid for temperatures over 80F: if T < 80.0 or R < 40.0: return T hi_F = -42.379 + 2.04901523 * T + 10.14333127 * R - 0.22475541 * T * R - 6.83783e-3 * T ** 2\ - 5.481717e-2 * R ** 2 + 1.22874e-3 * T ** 2 * R + 8.5282e-4 * T * R ** 2 - 1.99e-6 * T ** 2 * R ** 2 if hi_F < T: hi_F = T return hi_F def heatindexC(T_C, R): if T_C is None or R is None: return None T_F = CtoF(T_C) hi_F = heatindexF(T_F, R) return FtoC(hi_F) def heating_degrees(t, base): return max(base - t, 0) if t is not None else None def cooling_degrees(t, base): return max(t - base, 0) if t is not None else None def altimeter_pressure_US(SP_inHg, Z_foot, algorithm='aaASOS'): """Calculate the altimeter pressure, given the raw, station pressure in inHg and the altitude in feet. Examples: >>> print "%.2f" % altimeter_pressure_US(28.0, 0.0) 28.00 >>> print "%.2f" % altimeter_pressure_US(28.0, 1000.0) 29.04 """ if SP_inHg is None or Z_foot is None: return None if SP_inHg <= 0.008859: return None return weewx.uwxutils.TWxUtilsUS.StationToAltimeter(SP_inHg, Z_foot, algorithm=algorithm) def altimeter_pressure_Metric(SP_mbar, Z_meter, algorithm='aaASOS'): """Convert from (uncorrected) station pressure to altitude-corrected pressure. Examples: >>> print "%.1f" % altimeter_pressure_Metric(948.08, 0.0) 948.2 >>> print "%.1f" % altimeter_pressure_Metric(948.08, 304.8) 983.4 """ if SP_mbar is None or Z_meter is None: return None if SP_mbar <= 0.3: return None return weewx.uwxutils.TWxUtils.StationToAltimeter(SP_mbar, Z_meter, algorithm=algorithm) def _etterm(elev_meter, t_C): """Calculate elevation/temperature term for sea level calculation.""" t_K = CtoK(t_C) return math.exp(-elev_meter / (t_K * 29.263)) def sealevel_pressure_Metric(sp_mbar, elev_meter, t_C): """Convert station pressure to sea level pressure. This implementation was copied from wview. sp_mbar - station pressure in millibars elev_meter - station elevation in meters t_C - temperature in degrees Celsius bp - sea level pressure (barometer) in millibars """ if sp_mbar is None or elev_meter is None or t_C is None: return None pt = _etterm(elev_meter, t_C) bp_mbar = sp_mbar / pt if pt != 0 else 0 return bp_mbar def sealevel_pressure_US(sp_inHg, elev_foot, t_F): if sp_inHg is None or elev_foot is None or t_F is None: return None sp_mbar = sp_inHg / INHG_PER_MBAR elev_meter = elev_foot * METER_PER_FOOT t_C = FtoC(t_F) slp_mbar = sealevel_pressure_Metric(sp_mbar, elev_meter, t_C) slp_inHg = slp_mbar * INHG_PER_MBAR return slp_inHg def calculate_rain(newtotal, oldtotal): """Calculate the rain differential given two cumulative measurements.""" if newtotal is not None and oldtotal is not None: if newtotal >= oldtotal: delta = newtotal - oldtotal else: syslog.syslog(syslog.LOG_INFO, "wxformulas: rain counter reset detected: new=%s old=%s" % (newtotal, oldtotal)) delta = None else: delta = None return delta def solar_rad_Bras(lat, lon, altitude_m, ts=None, nfac=2): """Calculate maximum solar radiation using Bras method http://www.ecy.wa.gov/programs/eap/models.html lat, lon - latitude and longitude in decimal degrees altitude_m - altitude in meters ts - timestamp as unix epoch nfac - atmospheric turbidity (2=clear, 4-5=smoggy) Example: >>> for t in range(0,24): ... print "%.2f" % solar_rad_Bras(42, -72, 0, t*3600+1422936471) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.86 100.81 248.71 374.68 454.90 478.76 443.47 353.23 220.51 73.71 0.00 0.00 0.00 0.00 0.00 0.00 """ from weewx.almanac import Almanac if ts is None: ts = time.time() sr = 0.0 try: alm = Almanac(ts, lat, lon, altitude_m) el = alm.sun.alt # solar elevation degrees from horizon R = alm.sun.earth_distance # NREL solar constant W/m^2 nrel = 1367.0 # radiation on horizontal surface at top of atmosphere (bras eqn 2.9) sinel = math.sin(math.radians(el)) io = sinel * nrel / (R * R) if sinel >= 0: # optical air mass (bras eqn 2.22) m = 1.0 / (sinel + 0.15 * math.pow(el + 3.885, -1.253)) # molecular scattering coefficient (bras eqn 2.26) a1 = 0.128 - 0.054 * math.log(m) / math.log(10.0) # clear-sky radiation at earth surface W / m^2 (bras eqn 2.25) sr = io * math.exp(-nfac * a1 * m) except (AttributeError, ValueError, OverflowError): sr = None return sr def solar_rad_RS(lat, lon, altitude_m, ts=None, atc=0.8): """Calculate maximum solar radiation Ryan-Stolzenbach, MIT 1972 http://www.ecy.wa.gov/programs/eap/models.html lat, lon - latitude and longitude in decimal degrees altitude_m - altitude in meters ts - time as unix epoch atc - atmospheric transmission coefficient (0.7-0.91) Example: >>> for t in range(0,24): ... print "%.2f" % solar_rad_RS(42, -72, 0, t*3600+1422936471) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.09 79.31 234.77 369.80 455.66 481.15 443.44 346.81 204.64 52.63 0.00 0.00 0.00 0.00 0.00 0.00 """ from weewx.almanac import Almanac if atc < 0.7 or atc > 0.91: atc = 0.8 if ts is None: ts = time.time() sr = 0.0 try: alm = Almanac(ts, lat, lon, altitude_m) el = alm.sun.alt # solar elevation degrees from horizon R = alm.sun.earth_distance z = altitude_m nrel = 1367.0 # NREL solar constant, W/m^2 sinal = math.sin(math.radians(el)) if sinal >= 0: # sun must be above horizon rm = math.pow((288.0 - 0.0065 * z) / 288.0, 5.256) / (sinal + 0.15 * math.pow(el + 3.885, -1.253)) toa = nrel * sinal / (R * R) sr = toa * math.pow(atc, rm) except (AttributeError, ValueError, OverflowError): sr = None return sr def wetbulp(t_C, rh): """Calculate the wet bulp temperatur t_C - temperature in degrees Celsius rh - relative humidity [0-100] """ dp_ce = dewpointC(t_C, rh) if dp_ce is None: return None #wbt2 = t_C * math.atan(0.151977 * math.pow((rh + 8.313659), 0.5)) + math.atan(t_C + rh) - math.atan(rh - 1.676331) + 0.00391838 * math.pow(rh, 1.5) * math.atan(0.023101 * rh) - 4.686035 wbt = -5.809 + 0.058 * rh + 0.697 * t_C + 0.003 * rh * t_C return wbt def cloudbase_Metric(t_C, rh, altitude_m): """Calculate the cloud base in meters t_C - temperature in degrees Celsius rh - relative humidity [0-100] altitude_m - altitude in meters """ dp_C = dewpointC(t_C, rh) if dp_C is None: return None cb = (t_C - dp_C) * 1000 / 2.5 return altitude_m + cb * METER_PER_FOOT if cb is not None else None def cloudbase_US(t_F, rh, altitude_ft): """Calculate the cloud base in feet t_F - temperature in degrees Fahrenheit rh - relative humidity [0-100] altitude_ft - altitude in feet """ dp_F = dewpointF(t_F, rh) if dp_F is None: return None cb = altitude_ft + (t_F - dp_F) * 1000.0 / 4.4 return cb def humidexC(t_C, rh): """Calculate the humidex Reference (look under heading "Humidex"): http://climate.weather.gc.ca/climate_normals/normals_documentation_e.html?docID=1981 t_C - temperature in degree Celsius rh - relative humidity [0-100] Examples: >>> print "%.2f" % humidexC(30.0, 80.0) 43.64 >>> print "%.2f" % humidexC(30.0, 20.0) 30.00 >>> print "%.2f" % humidexC(0, 80.0) 0.00 >>> print humidexC(30.0, None) None """ try: dp_C = dewpointC(t_C, rh) dp_K = CtoK(dp_C) e = 6.11 * math.exp(5417.7530 * (1 / 273.16 - 1 / dp_K)) h = 0.5555 * (e - 10.0) except (ValueError, OverflowError, TypeError): return None return t_C + h if h > 0 else t_C def humidexF(t_F, rh): """Calculate the humidex in degree Fahrenheit t_F - temperature in degree Fahrenheit rh - relative humidity [0-100] """ if t_F is None: return None h_C = humidexC(FtoC(t_F), rh) return CtoF(h_C) if h_C is not None else None def apptempC(t_C, rh, ws_mps): """Calculate the apparent temperature in degree Celsius t_C - temperature in degree Celsius rh - relative humidity [0-100] ws_mps - wind speed in meters per second http://www.bom.gov.au/info/thermal_stress/#atapproximation AT = Ta + 0.33*e - 0.70*ws - 4.00 where AT and Ta (air temperature) are deg-C e is water vapor pressure ws is wind speed (m/s) at elevation of 10 meters e = rh / 100 * 6.105 * exp(17.27 * Ta / (237.7 + Ta)) rh is relative humidity http://www.ncdc.noaa.gov/societal-impacts/apparent-temp/ AT = -2.7 + 1.04*T + 2.0*e -0.65*v where AT and T (air temperature) are deg-C e is vapor pressure in kPa v is 10m wind speed in m/sec """ if t_C is None: return None if rh is None or rh < 0 or rh > 100: return None if ws_mps is None or ws_mps < 0: return None try: e = (rh / 100.0) * 6.105 * math.exp(17.27 * t_C / (237.7 + t_C)) at_C = t_C + 0.33 * e - 0.7 * ws_mps - 4.0 except (ValueError, OverflowError): at_C = None return at_C def apptempF(t_F, rh, ws_mph): """Calculate apparent temperature in degree Fahrenheit t_F - temperature in degree Fahrenheit rh - relative humidity [0-100] ws_mph - wind speed in miles per hour """ if t_F is None: return None if rh is None or rh < 0 or rh > 100: return None if ws_mph is None or ws_mph < 0: return None t_C = FtoC(t_F) ws_mps = ws_mph * METER_PER_MILE / 3600.0 at_C = apptempC(t_C, rh, ws_mps) return CtoF(at_C) if at_C is not None else None def beaufort(ws_kts): """Return the beaufort number given a wind speed in knots""" if ws_kts is None: return None elif ws_kts < 1: return 0 elif ws_kts < 4: return 1 elif ws_kts < 7: return 2 elif ws_kts < 11: return 3 elif ws_kts < 17: return 4 elif ws_kts < 22: return 5 elif ws_kts < 28: return 6 elif ws_kts < 34: return 7 elif ws_kts < 41: return 8 elif ws_kts < 48: return 9 elif ws_kts < 56: return 10 elif ws_kts < 64: return 11 return 12 def equation_of_time(doy): """Equation of time in minutes. Plus means sun leads local time. Example (1 October): >>> print "%.4f" % equation_of_time(274) 0.1889 """ b = 2 * math.pi * (doy - 81) / 364.0 return 0.1645 * math.sin(2 * b) - 0.1255 * math.cos(b) - 0.025 * math.sin(b) def hour_angle(t_utc, longitude, doy): """Solar hour angle at a given time in radians. t_utc: The time in UTC. longitude: the longitude in degrees doy: The day of year Returns hour angle in radians. 0 <= omega < 2*pi Example: >>> print "%.4f radians" % hour_angle(15.5, -16.25, 274) 0.6821 radians >>> print "%.4f radians" % hour_angle(0, -16.25, 274) 2.9074 radians """ Sc = equation_of_time(doy) omega = (math.pi / 12.0) * (t_utc + longitude / 15.0 + Sc - 12) if omega < 0: omega += 2.0 * math.pi return omega def solar_declination(doy): """Solar declination for the day of the year in radians Example (1 October is the 274th day of the year): >>> print "%.6f" % solar_declination(274) -0.075274 """ return 0.409 * math.sin(2.0 * math.pi * doy / 365 - 1.39) def sun_radiation(doy, latitude_deg, longitude_deg, tod_utc, interval): """Extraterrestrial radiation. Radiation at the top of the atmosphere doy: Day-of-year latitude_deg, longitude_deg: Lat and lon in degrees tod_utc: Time-of-day (UTC) at the end of the interval in hours (0-24) interval: The time interval over which the radiation is to be calculated in hours Returns the (average?) solar radiation over the time interval in MJ/m^2/hr Example: >>> print "%.3f" % sun_radiation(doy=274, latitude_deg=16.217, ... longitude_deg=-16.25, tod_utc=16.0, interval=1.0) 3.543 """ # Solar constant in MJ/m^2/hr Gsc = 4.92 delta = solar_declination(doy) earth_distance = 1.0 + 0.033 * math.cos(2.0 * math.pi * doy / 365.0) # dr start_utc = tod_utc - interval stop_utc = tod_utc start_omega = hour_angle(start_utc, longitude_deg, doy) stop_omega = hour_angle(stop_utc, longitude_deg, doy) latitude_radians = math.radians(latitude_deg) part1 = (stop_omega - start_omega) * math.sin(latitude_radians) * math.sin(delta) part2 = math.cos(latitude_radians) * math.cos(delta) * (math.sin(stop_omega) - math.sin(start_omega)) # http://www.fao.org/docrep/x0490e/x0490e00.htm Eqn 28 Ra = (12.0 / math.pi) * Gsc * earth_distance * (part1 + part2) if Ra < 0: Ra = 0 return Ra def longwave_radiation(Tmin_C, Tmax_C, ea, Rs, Rso, rh): """Calculate the net long-wave radiation. Ref: http://www.fao.org/docrep/x0490e/x0490e00.htm Eqn 39 Tmin_C: Minimum temperature during the calculation period Tmax_C: Maximum temperature during the calculation period ea: Actual vapor pressure in kPa Rs: Measured radiation. See below for units. Rso: Calculated clear-wky radiation. See below for units. rh: Relative humidity in percent Because the formula uses the ratio of Rs to Rso, their actual units do not matter, so long as they use the same units. Returns back radiation in MJ/m^2/day Example: >>> print "%.1f mm/day" % longwave_radiation(Tmin_C=19.1, Tmax_C=25.1, ea=2.1, Rs=14.5, Rso=18.8, rh=50) 3.5 mm/day Night time example. Set rh = 40% to reproduce the Rs/Rso ratio of 0.8 used in the paper. >>> print "%.1f mm/day" % longwave_radiation(Tmin_C=28, Tmax_C=28, ea=3.402, Rs=0, Rso=0, rh=40) 2.4 mm/day """ # Calculate temperatures in Kelvin: Tmin_K = Tmin_C + 273.16 Tmax_K = Tmax_C + 273.16 # Stefan-Boltzman constant in MJ/K^4/m^2/day sigma = 4.903e-09 # Use the ratio of measured to expected radiation as a measure of cloudiness, but # only if it's daylight if Rso: cloud_factor = Rs / Rso else: # If it's nighttime (no expected radiation), then use this totally made up formula if rh > 80: # Humid. Lots of clouds cloud_factor = 0.3 elif rh > 40: # Somewhat humid. Modest cloud cover cloud_factor = 0.5 else: # Low humidity. No clouds. cloud_factor = 0.8 # Calculate the longwave (back) radiation (Eqn 39). Result will be in MJ/m^2/day. Rnl_part1 = sigma * (Tmin_K ** 4 + Tmax_K ** 4) / 2.0 Rnl_part2 = (0.34 - 0.14 * math.sqrt(ea)) Rnl_part3 = (1.35 * cloud_factor - 0.35) Rnl = Rnl_part1 * Rnl_part2 * Rnl_part3 return Rnl def evapotranspiration_Metric(Tmin_C, Tmax_C, rh_min, rh_max, sr_mean_wpm2, ws_mps, wind_height_m, latitude_deg, longitude_deg, altitude_m, timestamp): """Calculate the rate of evapotranspiration during a one hour time period. Ref: http://www.fao.org/docrep/x0490e/x0490e00.htm. (The document http://edis.ifas.ufl.edu/ae459 is also helpful) Tmin_C: Minimum temperature during the hour in degrees Celsius Tmax_C: Maximum temperature during the hour in degrees Celsius rh_min: Minimum relative humidity during the hour in percent. rh_max: Maximum relative humidity during the hour in percent. sr_mean_wpm2: Mean solar radiation during the hour in watts per sq meter ws_mps: Average wind speed during the hour in meters per second wind_height_m: Height in meters at which windspeed is measured latitude_deg, longitude_deg: Latitude, longitude of the station in degrees altitude_m: Altitude of the station in meters. timestamp: The time, as unix epoch time, at the end of the hour. Returns: Evapotranspiration in mm/hr Example (Example 19 in the reference document): >>> sr_mean_wpm2 = 680.56 # == 2.45 MJ/m^2/hr >>> timestamp = 1475337600 # 1-Oct-2016 at 16:00UTC >>> print "ET0 = %.2f mm/hr" % evapotranspiration_Metric(Tmin_C=38, Tmax_C=38, rh_min=52, rh_max=52, ... sr_mean_wpm2=sr_mean_wpm2, ws_mps=3.3, wind_height_m=2, ... latitude_deg=16.217, longitude_deg=-16.25, altitude_m=8, timestamp=timestamp) ET0 = 0.63 mm/hr Another example, this time for night >>> sr_mean_wpm2 = 0.0 # night time >>> timestamp = 1475294400 # 1-Oct-2016 at 04:00UTC (0300 local) >>> print "ET0 = %.2f mm/hr" % evapotranspiration_Metric(Tmin_C=28, Tmax_C=28, rh_min=90, rh_max=90, ... sr_mean_wpm2=sr_mean_wpm2, ws_mps=3.3, wind_height_m=2, ... latitude_deg=16.217, longitude_deg=-16.25, altitude_m=8, timestamp=timestamp) ET0 = 0.03 mm/hr """ if None in (Tmin_C, Tmax_C, rh_min, rh_max, sr_mean_wpm2, ws_mps, latitude_deg, longitude_deg, timestamp): return None if wind_height_m is None: wind_height_m = 2.0 if altitude_m is None: altitude_m = 0.0 # Numerator and denominator terms for the reference crop type cn = 37 cd = 0.34 # Albedo. for grass reference crop albedo = 0.23 # figure out the day of year [1-366] from the timestamp doy = time.localtime(timestamp)[7] - 1 # Calculate the UTC time-of-day in hours time_tt_utc = time.gmtime(timestamp) tod_utc = time_tt_utc.tm_hour + time_tt_utc.tm_min / 60.0 + time_tt_utc.tm_sec / 3600.0 # Calculate mean temperature tavg_C = (Tmax_C + Tmin_C) / 2.0 # Mean humidity rh_avg = (rh_min + rh_max) / 2.0 # Adjust windspeed for height u2 = 4.87 * ws_mps / math.log(67.8 * wind_height_m - 5.42) # Calculate the atmospheric pressure in kPa p = 101.3 * math.pow((293.0 - 0.0065 * altitude_m) / 293.0, 5.26) # Calculate the psychrometric constant in kPa/C (Eqn 8) gamma = 0.665e-03 * p # Calculate mean saturation vapor pressure, converting from hPa to kPa (Eqn 12) etmin = weewx.uwxutils.TWxUtils.SaturationVaporPressure(Tmin_C, 'vaTeten') / 10.0 etmax = weewx.uwxutils.TWxUtils.SaturationVaporPressure(Tmax_C, 'vaTeten') / 10.0 e0T = (etmin + etmax) / 2.0 # Calculate the slope of the saturation vapor pressure curve in kPa/C (Eqn 13) delta = 4098.0 * (0.6108 * math.exp(17.27 * tavg_C / (tavg_C + 237.3))) / \ ((tavg_C + 237.3) * (tavg_C + 237.3)) # Calculate actual vapor pressure from relative humidity (Eqn 17) ea = (etmin * rh_max + etmax * rh_min) / 200.0 # Convert solar radiation from W/m^2 to MJ/m^2/hr Rs = sr_mean_wpm2 * 3.6e-3 # Net shortwave (measured) radiation in MJ/m^2/hr (eqn 38) Rns = (1.0 - albedo) * Rs # Extraterrestrial radiation in MJ/m^2/hr Ra = sun_radiation(doy, latitude_deg, longitude_deg, tod_utc, interval=1.0) # Clear sky solar radiation in MJ/m^2/hr (eqn 37) Rso = (0.75 + 2e-5 * altitude_m) * Ra # Longwave (back) radiation. Convert from MJ/m^2/day to MJ/m^2/hr (Eqn 39): Rnl = longwave_radiation(Tmin_C, Tmax_C, ea, Rs, Rso, rh_avg) / 24.0 # Calculate net radiation at the surface in MJ/m^2/hr (Eqn. 40) Rn = Rns - Rnl # Calculate the soil heat flux. (see section "For hourly or shorter # periods" in http://www.fao.org/docrep/x0490e/x0490e07.htm#radiation G = 0.1 * Rn if Rs else 0.5 * Rn # Put it all together. Result is in mm/hr (Eqn 53) ET0 = (0.408 * delta * (Rn - G) + gamma * (cn / (tavg_C + 273)) * u2 * (e0T - ea)) / (delta + gamma * (1 + cd * u2)) # We don't allow negative ET's if ET0 < 0: ET0 = 0 return ET0 def evapotranspiration_US(Tmin_F, Tmax_F, rh_min, rh_max, sr_mean_wpm2, ws_mph, wind_height_ft, latitude_deg, longitude_deg, altitude_ft, timestamp): """Calculate the rate of evapotranspiration during a one hour time period, returning result in inches/hr. Tmin_F: Minimum temperature during the hour in degrees Fahrenheit Tmax_F: Maximum temperature during the hour in degrees Fahrenheit rh_min: Minimum relative humidity during the hour in percent. rh_max: Maximum relative humidity during the hour in percent. sr_mean_wpm2: Mean solar radiation during the hour in watts per sq meter ws_mph: Average wind speed during the hour in miles per hour wind_height_ft: Height in feet at which windspeed is measured latitude_deg, longitude_deg: Latitude, longitude of the station in degrees altitude_ft: Altitude of the station in feet. timestamp: The time, as unix epoch time, at the end of the hour. Returns: Evapotranspiration in inches/hr Example (using data from HR station): >>> sr_mean_wpm2 = 860 >>> timestamp = 1469829600 # 29-July-2016 22:00 UTC (15:00 local time) >>> print "ET0 = %.3f in/hr" % evapotranspiration_US(Tmin_F=87.8, Tmax_F=89.1, rh_min=34, rh_max=38, ... sr_mean_wpm2=sr_mean_wpm2, ws_mph=9.58, wind_height_ft=6, ... latitude_deg=45.7, longitude_deg=-121.5, altitude_ft=700, timestamp=timestamp) ET0 = 0.028 in/hr """ try: Tmin_C = FtoC(Tmin_F) Tmax_C = FtoC(Tmax_F) ws_mps = ws_mph * METER_PER_MILE / 3600.0 wind_height_m = wind_height_ft * METER_PER_FOOT altitude_m = altitude_ft * METER_PER_FOOT except TypeError: return None evt = evapotranspiration_Metric(Tmin_C=Tmin_C, Tmax_C=Tmax_C, rh_min=rh_min, rh_max=rh_max, sr_mean_wpm2=sr_mean_wpm2, ws_mps=ws_mps, wind_height_m=wind_height_m, latitude_deg=latitude_deg, longitude_deg=longitude_deg, altitude_m=altitude_m, timestamp=timestamp) return evt / MM_PER_INCH if evt is not None else None if __name__ == "__main__": import doctest if not doctest.testmod().failed: print("PASSED")
