From b4fa0f335064502cf5c92a75a9dcebcb555c134e Mon Sep 17 00:00:00 2001 From: Samer Albahra Date: Thu, 23 Apr 2015 14:38:14 -0500 Subject: [PATCH] Add support for calculating Evapotranspiration --- application.py | 137 +++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 126 insertions(+), 11 deletions(-) diff --git a/application.py b/application.py index 0a92dda..a10b88a 100644 --- a/application.py +++ b/application.py @@ -1,5 +1,5 @@ #!/usr/bin/python -import urllib, urllib2, cgi, re +import urllib, urllib2, cgi, re, math import json, datetime, time, sys, calendar import pytz, ephem from datetime import datetime, timedelta, date @@ -33,6 +33,18 @@ def isFloat(s): return 0 return 1 +def F2C(temp): + return (temp-32)*5/9 + +def C2F(temp): + return temp*9/5 + 32 + +def mm2in(x): + return x*0.03937008 + +def ft2m(x): + return x*0.3048 + def IP2Int(ip): o = map(int, ip.split('.')) res = (16777216 * o[0]) + (65536 * o[1]) + (256 * o[2]) + o[3] @@ -44,6 +56,78 @@ def getClientAddress(environ): except KeyError: return environ['REMOTE_ADDR'] +def computeETs(latitude, longitude, elevation, temp_high, temp_low, temp_avg, hum_high, hum_low, hum_avg, wind, solar): + tm = time.gmtime() + dayofyear = tm.tm_yday + + latitude = safe_float(latitude, 0); + longitude = safe_float(longitude, 0); + + # Converted values + El = ft2m(elevation) + Rs = float(solar) * 0.0864 # W/m2 to MJ/d /m2 + Tx = F2C(float(temp_high)) + Tn = F2C(float(temp_low)) + Tm = F2C(float(temp_avg)) + RHx = float(hum_high) + RHn = float(hum_low) + RHm = float(hum_avg) + Td = Tm-(100-RHm)/5 # approx. dewpoint (daily mean) + U2 = float(wind) * 0.44704 # wind speed in m/s + + # Step 1: Extraterrestrial radiation + + Gsc = 0.082 + sigma = 4.90e-9 + phi = math.pi * latitude / 180 + dr = 1 + 0.033*math.cos(2*math.pi*dayofyear/365) + delta = 0.409*math.sin(2*math.pi*dayofyear/365 - 1.39) + omegas = math.acos(-math.tan(phi)*math.tan(delta)) + Ra = (24*60/math.pi)*Gsc*dr*(omegas*math.sin(delta)*math.sin(phi)+math.cos(phi)*math.cos(delta)*math.sin(omegas)) + + # Step 2: Daily net radiation + + Rso = Ra*(0.75 + 2.0e-5 * El) # 5 + Rns = (1-0.23)*Rs + f = 1.35*Rs/Rso - 0.35 # 7 + + esTx = 0.6108*math.exp(17.27*Tx/(Tx+237.3)) # 8 + esTn = 0.6108*math.exp(17.27*Tn/(Tn+237.3)) + ed = (esTx*RHn/100 + esTn*RHx/100)/2 # 10 + ea = (esTx+esTn)/2 # 22 + + epsilonp = 0.34-0.14*math.sqrt(ea) # 12 + Rnl = -f*epsilonp*sigma*((Tx+273.14)**4 + (Tn+273.15)**4)/2 # 13 + Rn = Rns + Rnl + + # Step 3: variables needed to compute ET + + beta = 101.3*((293-0.0065*El)/293)**5.26 # 15 + lam = 2.45 + gamma = 0.00163*beta/lam + e0 = 0.6108*math.exp(17.27*Tm/(Tm+237.3)) # 19 + Delta = 4099*e0/(Tm+237.3)**2 # 20 + G = 0 + ea = (esTx+esTn)/2 + + # Step 4: calculate ETh + + ETh = 0.408*(0.0023*Ra*(Tm+17.8)*math.sqrt(Tx-Tn)) # 23 + + # Step 5: calculate ET0 + + R0 = 0.408*Delta*(Rn-G)/(Delta+gamma*(1+0.34*U2)) # 24 + A0 = (900*gamma/(Tm+273))*U2*(ea-ed) / (Delta+gamma*(1+0.34*U2)) # 25 + ET0 = R0 + A0 + + # Step 6: calculate ETr + + Rr = 0.408*Delta*(Rn-G) / (Delta+gamma*(1+0.38*U2)) # 27 + Ar = (1600*gamma/(Tm+273))*U2*(ea-ed) / (Delta+gamma*(1+0.38*U2)) # 28 + ETr = Rr + Ar + + return (mm2in(ETh), mm2in(ET0), mm2in(ETr)) + def application(environ, start_response): path = environ.get('PATH_INFO') uwt = re.match('/weather(\d+)\.py',path) @@ -72,8 +156,8 @@ def application(environ, start_response): else: fwv = '' - restrict, maxh, minh, meant, pre, pre_today, h_today, sunrise, sunset, scale, toffset = [0, -1, -1, -500, -1, -1, -1, -1, -1, -1, -1] - + solar, wind, avehumidity, minhumidity, maxhumidity, maxt, mint, elevation, restrict, maxh, minh, meant, pre, pre_today, h_today, sunrise, sunset, scale, toffset = [0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -500, -1, -1, -1, -1, -1, -1, -1] + ET = [0,0,0] eip = IP2Int(getClientAddress(environ)) # if loc is GPS coordinate itself @@ -94,30 +178,55 @@ def application(environ, start_response): # if loc is pws, query wunderground geolookup to get GPS coordinates if loc.startswith('pws:') or loc.startswith('icao:'): try: - req = urllib2.urlopen('http://api.wunderground.com/api/'+key+'/geolookup/q/'+urllib.quote(loc)+'.json') + req = urllib2.urlopen('http://api.wunderground.com/api/'+key+'/conditions/forecast/q/'+urllib.quote(loc)+'.json') dat = json.load(req) - if dat['location']: - v = dat['location']['lat'] + if 'current_observation' in dat: + v = dat['current_observation']['observation_location']['latitude'] if v and isFloat(v): lat = v - v = dat['location']['lon'] + v = dat['current_observation']['observation_location']['longitude'] if v and isFloat(v): lon = v - v = dat['location']['tz_long'] + v = dat['current_observation']['observation_location']['elevation'] + if v: + elevation = safe_int(int(v.split()[0]), 0) + v = dat['current_observation']['solarradiation'] + if v: + solar = safe_int(v, 0) + v = dat['current_observation']['local_tz_long'] if v: tzone = v else: - v = dat['location']['tz'] + v = dat['current_observation']['local_tz_short'] if v: tzone = v + forecast = dat['forecast']['simpleforecast']['forecastday'][0] + + v = forecast['high']['fahrenheit'] + if v: + maxt = safe_int(v, 0) + v = forecast['low']['fahrenheit'] + if v: + mint = safe_int(v, 0) + v = forecast['avehumidity'] + if v: + avehumidity = safe_int(v, 0) + v = forecast['maxhumidity'] + if v: + maxhumidity = safe_int(v, 0) + v = forecast['minhumidity'] + if v: + minhumidity = safe_int(v, 0) + v = forecast['avewind']['mph'] + if v: + wind = safe_int(v, 0) + except: lat = None lon = None tzone = None - #loc = loc.replace(' ','_') - # now do autocomplete lookup to get GPS coordinates if lat==None or lon==None: try: @@ -195,6 +304,7 @@ def application(environ, start_response): if v: h_today = safe_float(v, h_today) + # Check which weather method is being used if ((uwt & ~(1 << 7)) == 1): # calculate water time scale, per https://github.com/rszimm/sprinklers_pi/blob/master/Weather.cpp hf = 0 @@ -217,6 +327,11 @@ def application(environ, start_response): if (scale>200): scale = 200 + elif ((uwt & ~(1 << 7)) == 2): + ET = computeETs(lat, lon, elevation, maxt, mint, meant, maxhumidity, minhumidity, avehumidity, wind, solar) + # TODO: Actually generate correct scale using ET (ET[1] is ET0 for short canopy) + scale = safe_int(ET[1] * 100, -1) + # Check weather modifier bits and apply scale modification if ((uwt>>7) & 1): # California modification to prevent watering when rain has occured within 48 hours