Created
September 30, 2014 10:44
-
-
Save anonymous/7f52c4d04f1a5658dace to your computer and use it in GitHub Desktop.
Envirotranspiration script for python and Weatherunderground
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/python | |
import urllib, json, sys, math, time, os, smtplib; | |
from datetime import date, timedelta, datetime | |
from email.mime.text import MIMEText | |
# Location constant | |
GMTOFFSET = 5 | |
LATITUDE = 45.4214 | |
LONGITUDE = -75.6919 | |
COUNTRY = "Canada" | |
CITY = "Ottawa" | |
KEY = #Need to add your key here | |
#debug level. 0 gives a summary | |
level = 0 | |
#number of days to go back into the history to compute deficit | |
window = 21 | |
# Mapping of conditions to a level of cloud cover. These can be adjusted | |
# since they are all made up anyway | |
conditions = { | |
"Blowing Snow" :8, | |
"Clear" :0, | |
"Fog" :5, | |
"Haze" :2, | |
"Heavy Blowing Snow" :9, | |
"Heavy Fog" :9, | |
"Heavy Low Drifting Snow" :10, | |
"Heavy Rain" :10, | |
"Heavy Rain Showers" :10, | |
"Heavy Thunderstorms and Rain" :10, | |
"Light Drizzle" :10, | |
"Light Freezing Rain" :10, | |
"Light Ice Pellets" :10, | |
"Light Rain" :10, | |
"Light Rain Showers" :10, | |
"Light Snow" :10, | |
"Light Snow Grains" :10, | |
"Light Snow Showers" :10, | |
"Light Thunderstorms and Rain" :10, | |
"Low Drifting Snow" :10, | |
"Mist" :3, | |
"Mostly Cloudy" :8, | |
"Overcast" :10, | |
"Partial Fog" :2, | |
"Partly Cloudy" :5, | |
"Patches of Fog" :2, | |
"Rain" :10, | |
"Rain Showers" :10, | |
"Scattered Clouds" :4, | |
"Shallow Fog" :3, | |
"Snow" :10, | |
"Snow Showers" :10, | |
"Thunderstorm" :10, | |
"Thunderstorms and Rain" :10, | |
"Unknown" :5, | |
} | |
# Print an attribute | |
def printAttr(indata, name, uiname): | |
print " " + uiname + " " + str(indata[name]) | |
# Get forecast data for the city. Unused right now. | |
# Results could be used to suppress watering based on forecast rainfall | |
def getForecastData(): | |
forecastURL = 'http://api.wunderground.com/api/' + KEY + '/forecast/q/' + COUNTRY + '/' + CITY + '.json' | |
response = urllib.urlopen(forecastURL).read(); | |
data = json.loads(response) | |
print CITY + ', ' + COUNTRY | |
# Forecast day. Could make this a loop to lookahead even more | |
day = 0 | |
if (level > 3): | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day]['date'], "pretty", "Date") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day]['high'], "celsius", "High temp") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day]['low'], "celsius", "Low temp") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day], "maxhumidity", "High humidity") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day], "avehumidity", "Average humidity") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day], "minhumidity", "Low humidity") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day]['maxwind'], "kph", "Max wind") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day]['avewind'], "kph", "Average wind") | |
printAttr(data['forecast']['simpleforecast']['forecastday'][day]['qpf_allday'], "mm", "Daily rainfall") | |
forecast = float(data['forecast']['simpleforecast']['forecastday'][day]['qpf_allday']['mm']) | |
return forecast | |
# Returns a calculation of saturation vapour pressure based on temperature in degrees | |
def saturationVapourPressure(T): | |
return 0.6108 * math.exp((17.27 * T) / (T + 237.3)) | |
def getHistoricalData(forecast): | |
totalBalance = 0 | |
for day in range(window,-1,-1): | |
today = date.today() - timedelta(day) | |
datestring = today.strftime("%Y%m%d") | |
#response = urllib.urlopen(historyURL).read(); | |
try: | |
with open("water/" + datestring) as f: | |
for line in f: | |
irrigation = float(line) | |
if (irrigation == -1): | |
print " Reset balance" | |
totalBalance = 0 | |
else: | |
print " " + str(float(line)) + "mm of irrigation" | |
totalBalance += float(line) | |
except: | |
pass | |
if (day < 4): | |
try: | |
os.remove("data/" + datestring) | |
except OSError: | |
pass | |
try: | |
data = json.load(open("data/" + datestring)) | |
source = "file " + datestring | |
except: | |
historyURL = 'http://api.wunderground.com/api/'+ KEY + '/history_'+datestring+'/q/Canada/Ottawa.json' | |
time.sleep(10) | |
response = urllib.urlopen(historyURL).read() | |
cachefile = open("data/" + datestring, 'w') | |
cachefile.write(response) | |
cachefile.close() | |
data = json.loads(urllib.urlopen(historyURL).read()) | |
source = historyURL | |
thedate = date(int(data['history']['dailysummary'][0]['date']['year']), | |
int(data['history']['dailysummary'][0]['date']['mon']), | |
int(data['history']['dailysummary'][0]['date']['mday'])) | |
dayOfYear = thedate.timetuple().tm_yday | |
if (level > 0): | |
print 'Data for ' + CITY + ', ' + COUNTRY + ' from ' + source + " on " + str(thedate) | |
if (level > 2 ): | |
printAttr(data['history']['dailysummary'][0], "maxtempm", "High temp") | |
printAttr(data['history']['dailysummary'][0], "meantempm", "Average temp") | |
printAttr(data['history']['dailysummary'][0], "mintempm", "Low temp") | |
printAttr(data['history']['dailysummary'][0], "maxhumidity", "High humidity") | |
printAttr(data['history']['dailysummary'][0], "minhumidity", "Low humidity") | |
printAttr(data['history']['dailysummary'][0], "maxwspdm", "Max wind") | |
printAttr(data['history']['dailysummary'][0], "meanwindspdm", "Average wind") | |
printAttr(data['history']['dailysummary'][0], "minwspdm", "Min wind") | |
printAttr(data['history']['dailysummary'][0], "precipm", "Daily rainfall") | |
printAttr(data['history']['dailysummary'][0], "meanpressurem", "Average air pressure") | |
# Calculate solar radiation from location | |
totalSolarRadiation = 0 | |
totalClearSkyIsolation = 0 | |
sunnyHours = 0 | |
# Get the conditions for the 24 hours of the requested day | |
for hour in range(0,24): | |
# Sometimes data is missing for an hour. If we don't find data, cloud cover will stay at | |
# -1 | |
cloudCover = -1 | |
# Look through the historical data we have | |
for period in range(0, len(data['history']['observations'])): | |
# Look for our hour in the date | |
if (int(data['history']['observations'][period]['date']['hour']) == hour): | |
# If there are conditions in the data, get them, and find the percent cloud cover | |
# for this hour | |
if (data['history']['observations'][period]['conds']): | |
cloudCover = float(conditions[data['history']['observations'][period]['conds']])/10 | |
cloudCoverString = data['history']['observations'][period]['conds'] | |
break; | |
# If we didn't find any conditions for this hour, assume the same conditions as the | |
# previous hour. Sometimes we are missing early data, but this is usually at night | |
# anyway, so we are safe | |
if (cloudCover == -1): cloudCover = previousCloudCover | |
previousCloudCover = cloudCover | |
# If we have data | |
if (cloudCover != -1): | |
# Track the number of sunny hours | |
if (cloudCover < 0.5): | |
sunnyHours = sunnyHours + 1 | |
# Find out the angle of the sun was in the middle of this hour as a good | |
# estimate | |
gmtHour = hour + GMTOFFSET + 0.5 | |
fractionalDay = (360/365.25)*(dayOfYear+gmtHour/24) | |
f = math.radians(fractionalDay) | |
declination = 0.396372 - 22.91327 * math.cos(f) + 4.02543 * math.sin(f) - 0.387205 * math.cos(2 * f) + 0.051967 * math.sin(2 * f) - 0.154527 * math.cos(3 * f) + 0.084798 * math.sin(3 * f) | |
timeCorrection = 0.004297 + 0.107029 * math.cos(f) - 1.837877 * math.sin(f) - 0.837378 * math.cos(2*f) - 2.340475*math.sin(2*f) | |
solarHour = (gmtHour + 0.5 - 12)*15 + LONGITUDE + timeCorrection | |
if (solarHour < -180): solarHour = solarHour + 360 | |
if (solarHour > 180): solarHour = solarHour - 360 | |
solarFactor = math.sin(math.radians(LATITUDE))*math.sin(math.radians(declination))+math.cos(math.radians(LATITUDE))*math.cos(math.radians(declination))*math.cos(math.radians(solarHour)) | |
sunElevation = math.degrees(math.asin(solarFactor)) | |
clearSkyInsolation = 990 * math.sin(math.radians(sunElevation))-30 | |
if (clearSkyInsolation < 0): clearSkyInsolation = 0 | |
solarRadiation = clearSkyInsolation * (1 - 0.75*(math.pow(cloudCover,3.4))) | |
# Accumulate clear sky radiation and solar radiation on the ground | |
totalSolarRadiation += solarRadiation | |
totalClearSkyIsolation += clearSkyInsolation | |
# Convert from Wh / m^2 / d | |
radiationAtSurface = totalSolarRadiation * 3600 / 1000 / 1000 # MJ / m^2 / d | |
# m/s at 2m above ground | |
windspeed = float(data['history']['dailysummary'][0]['meanwindspdm']) * 1000 / 3600 * 0.748 | |
pressure = float(data['history']['dailysummary'][0]['meanpressurem']) / 10 # kPa | |
tempAvg = float(data['history']['dailysummary'][0]['meantempm']) # degrees C | |
tempMin = float(data['history']['dailysummary'][0]['mintempm']) # degrees C | |
tempMax = float(data['history']['dailysummary'][0]['maxtempm']) # degrees C | |
humidMax = float(data['history']['dailysummary'][0]['maxhumidity']) # degrees C | |
humidMin = float(data['history']['dailysummary'][0]['minhumidity']) # degrees C | |
if (level > 2): | |
print " ground windspeed " + str(windspeed) | |
print " Sunny hours " + str(sunnyHours) | |
D = 4098 * saturationVapourPressure(tempAvg) / math.pow(tempAvg + 237.3,2) | |
g = 0.665e-3 * pressure | |
es = (saturationVapourPressure(tempMin) + saturationVapourPressure(tempMax)) / 2 | |
ea = saturationVapourPressure(tempMin) * humidMax / 200 + saturationVapourPressure(tempMax) * humidMin / 200 | |
vaporPressDeficit = es - ea | |
if (level == 2): | |
print " Solar Radiation " + str(totalSolarRadiation) | |
print " Clear Sky Solar Radiation " + str(totalClearSkyIsolation) | |
if (level > 3): | |
print " D " + str(D) | |
print " g " + str(g) | |
print " es " + str(es) | |
print " ea " + str(ea) | |
print " Vapor Pressure Deficit " + str(vaporPressDeficit) | |
print " Radiation at surface " + str(radiationAtSurface) | |
ETo = ((0.408 * D * radiationAtSurface) + (g * 900 * windspeed * vaporPressDeficit) / (tempAvg + 273)) / (D + g * (1 + 0.34 * windspeed)) | |
rainfall = float(data['history']['dailysummary'][0]['precipm']) | |
dailyBalance = rainfall - ETo | |
totalBalance += dailyBalance | |
if (level > 0): | |
print " Enviro-Transpiration " + str(ETo) | |
print " Daily Balance " + str(dailyBalance) | |
print " " + str(day) + " day balance " + str(totalBalance) | |
if (level > 0): | |
print " Forecast " + str(forecast) | |
print " Balance over " + str(window) + " days = " + str(totalBalance + forecast) | |
totalBalance += forecast | |
if (totalBalance > -5): | |
try: | |
os.remove("need") | |
except: | |
pass | |
f = open('skip','w') | |
f.write(str(totalBalance)+'\n') | |
f.close() | |
else: | |
try: | |
os.remove("skip") | |
except: | |
pass | |
f = open('need','w') | |
f.write(str(totalBalance)+'\n') | |
f.close() | |
forecast = getForecastData() | |
getHistoricalData(forecast) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment