Add support for calculating Evapotranspiration

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']:
+            if 'current_observation' in dat:
-                v = dat['location']['lat']
+                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