Envirotranspiration

  • 3
  • Idea
  • Updated 2 years ago
I wrote this python script to calculate ETo for a given location and compute a watering deficit for irrigation system.

Any way to attach a file in this system?
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes

Posted 4 years ago

  • 3
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
#!/usr/bin/python

import urllib, json, sys, math, time;
from datetime import date, timedelta, datetime

# Location constant
GMTOFFSET = 5
LATITUDE = 45.4214
LONGITUDE = -75.6919
COUNTRY = "Canada"
CITY = "Ottawa"
KEY = ""

#debug level. 0 gives a summary
level = 0

#number of days to go back into the history to compute deficit
window = 180

# 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(url, location):
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

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

# 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():

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:
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

getHistoricalData()
#getForecastData()
Photo of projectgreen

projectgreen

  • 2 Posts
  • 2 Reply Likes
Very cool! Thank you!
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
great piece of work!

A couple questions: Have you ever found such a script written visual basic?

And would the python script have to run on a raspberry pi?
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
I have never found such a script in any language, which is why I wrote this one. Not sure how easy it would be to translate to visual basic. Python has a lot of support for web API handling so that's why I did it that way. I use this on a Linux server to trigger a watering script every morning to see how much watering is needed. Could work on a pi too.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Thanks!

I control my sprinklers with a Webcontrol board that is controlled by Homeseer 2.5 on Windows XP. I just installed a program called ActivePython-2.5.5.7. I've read that a python script can be run from Homeseer with ActivePython running in the background (maybe). So I thought I would try your script in that manner.

I did find a JavaScript ET script. But it was 'written for' Wenatchee, WA. The solar radiation function coding has factors for those coordinates hard coded - no easy way to convert for different coordinates. Someone more familiar with that coding and the algorithm used could figure out how to make the coordinates variables like you did. I will try to paste the link here:
http://www.tfrec.wsu.edu/orchard/pet/...
The source is viewable. It's unfortunate the author didn't allow the coordinates to be selectable.

Thanks for your time!
Photo of Camden lowrance

Camden lowrance

  • 2 Posts
  • 0 Reply Likes
It would be much more helpful if you posted this as pre formatted python code, such as a github gist or pastebin.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
I agree. I'll see what I can do. I was hoping you could simply add an attachment here, but github is a good idea.
Photo of Camden lowrance

Camden lowrance

  • 2 Posts
  • 0 Reply Likes
Don't even bother making a repo, just paste the text into a gist:https://gist.github.com
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
https://gist.github.com/anonymous/e17...

Hope this works and helps.

I have not yet tested these results for sanity against other ET scripts, so if the results are different or you find problems, feel free to let me know.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Almost have the script working in windows.

Having trouble around:

if (cloudCover == -1): cloudCover = previousCloudCover
previousCloudCover = cloudCover

cloudCover is being set to previousCloudCover, but previousCloudCover is not defined at that point as far as I can tell.

Should the if statement end with a colon rather than how it now looks? It seems you lost some formatting, such as indentations, when you had to paste the script.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Can you try the one on the github gist? It has all the formatting and indents are very important in python.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
That looks a lot better.

One question: The variable, "source" initially refers to "file " + datestring

Then it is changed and set to historyURL, which includes the dataset I think. Is that OK?

I had to make a few changes to get the script to work in windows, and I do have it working now. But I'll change to the better formatted version.

You went to a lot of hard work on this. Thanks.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
It will take some time to modify this version to work.

First thing I find is a new folder "water". It is empty. I will have to determine whether a file will be populated there if it doesn't already exist.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Your Key is in the github version.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
The newer version script ran just fine. I changed it to my key. I used level 0 to get the summary. I assume the results are in mm.

The only problem remaining is the statement: with open("water/" + datestring) as f:

I've read that the only version of python to work with Homeseer is 2.5. And 2.5 does not recognize with open. I don't know how to handle the 'f:' without the with.

Thank you for your support and for the improvements you've made to this script. If you are curious, here is what the summary printed:

Ottawa, Canada
21 day balance -5.44131956622
20 day balance -10.503884595
19 day balance -12.16958391
18 day balance -12.8816373367
17 day balance -16.4811935658
16 day balance -9.9647276681
15 day balance -11.7920125933
14 day balance -13.7563728958
13 day balance -17.4464732077
12 day balance -15.8362476524
11 day balance -19.581386534
10 day balance -23.496282472
9 day balance -23.2542099084
8 day balance -13.9510685291
7 day balance -16.5776087087
6 day balance -19.001431717
5 day balance -22.3740392065
4 day balance -26.8574257149
3 day balance -31.6599764103
2 day balance -36.8340846739
1 day balance -41.7089125299
0 day balance -45.5322615361

Forecast 0.0
Balance over 21 days = -45.5322615361

Once I get the "with open(water..." section working, I will do comparisons with known data in my area.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
https://gist.github.com/anonymous/7f5...

This has the key removed and pulled to the top.

One open is mainly there so that the data you get from the site is cached and reused the next day. This minimized API usage, but could be removed easily enough.

The other open is to read the watering amounts from a separate file that is published from my controller. If you have a way to keep track of watering for the previous 21 days (or whatever), they you can adjust the deficit manually after the ET is done.

Did you change the GPS coordinates but keep the city name the same?

Mine ran this morning and this is what I got.

Ottawa, Canada
21 day balance -5.06256502876
20 day balance -6.72826434379
19 day balance -7.44031777047
18 day balance -11.0398739996
17 day balance -4.52340810188
16 day balance -6.35069302709
15 day balance -8.31505332957
14 day balance -12.0051536415
13 day balance -10.3949280862
12 day balance -14.1400669678
5.0mm of irrigation
11 day balance -13.0549629058
5.0mm of irrigation
10 day balance -7.81289034214
9 day balance 1.49025103707
5.0mm of irrigation
8 day balance 3.86371085749
7 day balance 1.43988784919
5.0mm of irrigation
6 day balance 3.06728035971
5.0mm of irrigation
5 day balance 3.58389385134
5.0mm of irrigation
4 day balance 3.78134315596
5.0mm of irrigation
3 day balance 3.60723489228
5.0mm of irrigation
2 day balance 3.73240703635
1 day balance -0.0909419698769
0 day balance -0.972172248588

Forecast 0.0
Balance over 21 days = -0.972172248588

No problem! Hope it works and is useful. Let me know of improvements or bugs.

Dave
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Thank you Dave. I appreciate your extra efforts to help.

I had changed the coordinates and the names to mine, but I put them back to yours to be sure my results compared with yours.

The only error I get is from the "with open" statement. The "open" one is fine.

I haven't tracked how much I have watered. But your script gives me a reason to. So I would plan to use the watering history file. I will also use the forecasted information to look forward to 'chance of rain' and high temps for a day (or two?).

I plan to incorporate this into Homeseer today to figure out how it can use this data to automatically adjust the duration of watering time.

Thanks again Dave,

Frank
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Hi Dave. The script is working fine in Homeseer HS2. I have it updating a variable in that program that could be accessed from another script. I've asked the author of the sprinkler controller software whether he could adapt his scripts to accept a value from your script that would adjust the Zone durations of a given schedule.

Would you let me know the format of the data within the water/ files?

I made as few changes as possible to your script, but still make it work with Python 2.5 and Homeseer. Also, it still works from my Windows command prompt without errors, so I assume it would also work in your Unix environment.

I have current/daily readings of ET and solar radiation in various counties here in Kansas. If you would like to have the links, let me know. I figure I will make comparisons between that information and the results from your script. But it seems the values yielded from your script look proper.

Thanks again for your hard work in writing the script.

Frank
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Some days show 24 sunny hours. It appears that factor isn't used in any calculations.

Maybe if the clearSkyInsolation is positive and there is no cloudcover = sunny day?
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Hmmmm... that was a sanity check only but it doesn't take into account whether the sun is actuallly up. ;) should be called "clearSkyHours" I guess.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Hmmm.. I guess somebody got ahold of my key from my other post and now I am going over limits. If anybody is using my key by accident, can we please check to fix that, or I'll have to pay or get a new key.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Nevermind. It was easy to regenerate, which I just did. If you script starts to fail, you will know why...
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
I'm glad some others are interested. You wrote a great script!

One question Dave, would you give a link to me of the formula or process you used to create the ET calculation? I would like to keep it with or in the script.

Still working to make it work correctly with Homeseer.

Thanks Dave,

Frank
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
http://www.fao.org/docrep/x0490e/x049...

This is the main source of data fir the equations in it.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
*for
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
http://www.shodor.org/os411/courses/_...


And this is what I used for the solar radiation equations.
(Edited)