Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Envirotranspiration script for python and Weatherunderground
#!/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")
print
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
print
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
You can’t perform that action at this time.