Evapotranspiration data from a Vantage Pro weather station

I've read how Berkinet and others are tweaking their watering schedules using evapotranspiration data from web sites like CIMIS, and been intrigued. Sadly where I live in Houston weather is highly localized and variable, and I didn't have a reliable source of nearby evapotranspiration data.

Then I got a Davis Vantage Pro 2 weather station for Father's day last June. I noticed the console displayed evapotranspiration data, but WeatherSnoop (which is otherwise awesome) didn't support it. I started doing a little digging in to how hard it would be to calculate it myself, and found a gold mine. Turns out the American Society of Civil Engineers decided to standardize evapotranspiration calculations about 10 years ago, and while the final report costs $30 to download, the draft (and all the appendices!) are downloadable for free from the University of Idaho.

Given how variable our local weather is, I opted to calculate the reference evapotranspiration constant based on hourly averages of weather data. I've had the thing working for a few days now, and it seems to agree with console calcs and common sense. If anyone who has a weather station that measures temperature, relative humidity, wind speed, and solar radiation wants to give it a try here it is, with apologies for my sad lack of Mad Python Skillz.

Code: Select all

"""
Calculate hourly evapotranspiration data from a weather station equipped to measure:
Relative Humidity, Temperature, Wind Speed (e.g., Davis Vantage Pro 2)
Relies on a separate script to totalize  relevant data from  device "WeatherSnoop" at
1 minute intervals.

This script runs on the hour, calculates ETo (Evapotranspiration reference value)
for Use by irrigation scripts
The rational for conversion of watts per square meter to
and from megajoules per square meter per hour are covered here:
http://rredc.nrel.gov/solar/pubs/sbf/b.html
In brief: An average of 1 watt for a 1 hour period = 3.6 kJ/h

These calculations rely on the approach outlined by
THE ASCE STANDARDIZED REFERENCE EVAPOTRANSPIRATION EQUATION
developed by the American Society of Civil Engineers, specifically those
for hourly data, and references to equations (e.g., eq. 55) refer to the
equation numbers used in that document.
Draft versions of the ASCE report and the appendices are downloadable at no charge from:
http://www.kimberly.uidaho.edu/water/asceewri/
"""
import math
from math import exp
import time

# Global constants
PI = 3.14159265
alt = 11 #in meters
lat = 29.568 #in decimal degrees
long = 95.05 #in decimal degrees
LZ = 90 # set longitude for center of time zone
Zw = 2 # height of wind speed measurement in meters above ground level
GSC = 4.92 # Global Solar Constant in MJ m^-2 h^-1
sigma = 2.042e-10 # Stefan-Boltzmann constant (in MJ K^-4 m^-2 h^-1)
albedo = 0.23 # generally accepted value for canopy reflection coefficient
Cn = 37 # numerator constant that for grass/hourly time steps (Table 1)
Cd_day = 0.24
Cd_night = 0.96


# Load totals for the last hourly measurement period from Indigo variables
totalSolRad = float(indigo.variables[1514345125].value) # "totalSolRad"
totalTempC = float(indigo.variables[900917121].value) # "totalTempC"
totalHumid  = float(indigo.variables[137383963].value) # "totalHumid"
totalPress = float(indigo.variables[1583976].value) # "totalPress"
totalWind  = float(indigo.variables[27494067].value) # "totalWind"
totalSampleCount = int(indigo.variables[582955114].value) # "totalSampleCount"
dailyETo = float(indigo.variables[1849325147].value) # "dailyETo"
runningETo = float(indigo.variables[1374037415].value) # "runningETo"

# Read historical cloud factor in case sun is too low in sky to accurately
# calculate a new one from solar radiation data
Fcd = float(indigo.variables[76116845].value) # "cloudFactor"

#  Calculate average values for input variables
avgSolRad = totalSolRad/totalSampleCount
avgTempC = totalTempC/totalSampleCount
avgHumid = totalHumid/totalSampleCount
avgPress = totalPress/totalSampleCount
avgWind = totalWind/totalSampleCount
indigo.server.log("avgSolRad = " + str(avgSolRad) + "W/m2")
indigo.server.log("avgTempC = " + str(avgTempC))
indigo.server.log("avgHumid = " + str(avgHumid))
indigo.server.log("avgPress = " + str(avgPress))
indigo.server.log("avgWind = " + str(avgWind))


# Calculate misc. conversions and date/time related info we will need
avgWind = 0.447 * avgWind # convert wind speed from MPH to m/s
avgTempK = avgTempC + 273.15 # Convert temperature to Kelvins
lt = time.localtime(time.time())
doy = lt[7] #Day of Year from time tuple
Rs = avgSolRad * 3.6e-3 # convert average wattage to megajoules per square meter per hour
indigo.server.log("Rs = " + str(Rs) + "MJ/m2h-1")

avgTempK = avgTempC + 273.15
P = 101.3 - 0.0115 * alt +5.44 * 10e-7 * alt**2 #atmospheric pressure from eq. 34
indigo.server.log("pressure = " + str(P))

# Correct wind speed measurement height to 2 meters using eq. 67
U2 = avgWind * (4.87/math.log(67.8 * Zw -5.42))

gamma = 0.000665 * P #Psychrometric constant from eq. 35
es = 0.6108 * exp(avgTempC*17.27/(avgTempC+237.3)) # saturation vapor pressure eq. 37
indigo.server.log("Psychrometric constant (gamma) = " + str(gamma))

# Using relative humidity data to calculate ea since weather station measures this
# and it is a "preferred method" as per ASCE
ea = es * avgHumid/100 # actual vapor pressure from eq. 41

# Calculate slope of saturation vapor pressure curve (delta) using eq. 36
delta = 2503 * exp(17.27 * avgTempC / (avgTempC + 237.3))/((avgTempC + 237.3)**2)
indigo.server.log("slope of vapor pressure (delta) = " + str(delta))

"""
this part gets a little ugly, because we have to figure out where the sun is in the sky
during the observation period, and deal with cases where either sunrise or sunset occur
during the observation period. We also need to correct for seasonal variation in when "solar
noon" occurs, and for daylight savings time.
"""
sd = 0.409 * math.sin(((2 * PI / 365) * doy - 1.39)) # solar declination from eq. 24
lat_rad = lat  * (PI / 180) # convert latitude from degrees to radians
sha = math.acos(-math.tan(lat_rad) * math.tan(sd)) # sunset hour angle from eq. 59
indigo.server.log("solar declination = " + str(sd))
indigo.server.log("sunset hour angle = " + str(sha))

# Calculate solar time angle at midpoint of time period
# correct for daylight savings time
t = lt[3] - 1  + lt[8] #Need to check if this works when NOT in Daylight Savings Time
t = t - 0.5 # use midpoint of measurement period

# Calculate seasonal corrections using eq. 57 & 58
b = 2 * PI * (doy - 81)/364
SC = 0.1645 * math.sin(2 * b) - 0.1255 * math.cos(b) - 0.025 * math.sin (b)

# Calculate sun angles for beginning and end of period
# Also deal with cases where sun rose or set during observation period
omega = PI/12 * ((t + 0.06667 * (LZ - long) + SC) - 12) #eq. 55
omega1 = omega - PI/24 #eq. 53
omega2 = omega + PI/24 #eq. 54
# Correct for sunset and sunrise using eq. 56
# recall that sunrise hour angle = -1 * sunset hour angle
if omega1 < -1 * sha:
   omega1 = -1 * sha
   omega1Dark = True
elif omega1 > sha:
   omega1 = sha
   omega1Dark = True
else:
   omega1Dark = False
   
if omega2 < -1 * sha:
   omega2 = -1 * sha
   omega2Dark = True
elif omega2 > sha:
   omega2 = sha
   omega2Dark = True
else:
   omega2Dark = False
   
if omega1 > omega2:
   omega1 = omega2
   
if omega1Dark:
   if omega2Dark:
      daytime = False
   else:
      daytime = True
else:
   daytime = True   

indigo.server.log("midpoint time t = " + str(t) + " hours")
indigo.server.log("seasonal correction = " + str(SC) + " hours")
indigo.server.log("omega = " + str(omega))
indigo.server.log("omega1 = " + str(omega1))
indigo.server.log("omega2 = " + str(omega2))
indigo.server.log("corrected sunset hour angle = " + str(sha))
indigo.server.log("daytime = " + str(daytime))

# Calculate sun angle above horizon
beta = math.asin(math.sin(lat_rad) * math.sin(sd) + math.cos(lat_rad) * math.cos(sd) * math.cos(omega))
indigo.server.log("current sun angle = " + str(beta) + " radians")
beta_deg = beta  * (180 / PI)
indigo.server.log("current angle = " + str(beta_deg) + " degrees")

# Calculate extraterrestrial radiation for period using eq. 48
# But only if some observations occurred during daylight hours
if daytime:
   IRD = 1 + (0.033 * math.cos((2 * PI / 365)* doy)) # inverse relative distance (earth to sun)
   Ra = (12/PI) * GSC * IRD * ((omega2 - omega1) * math.sin(lat_rad) * math.sin(sd) + math.cos(lat_rad) * math.cos(sd) * (math.sin(omega2) - math.sin (omega1)))
   indigo.server.log("Ra = " + str(Ra))

   #Calculate clear sky radiation using eq. 47
   Rso = Ra * (0.75 + 2e-5 * alt)
   indigo.server.log("Rso = " + str(Rso))

# Calculate cloudiness function using eq. 45
# But only if the sun is high enough in the sky to allow meaningful solar radiation measurement
if beta > 0.3:
   Fcd = (1.35 * (Rs/Rso)) - 0.35
   indigo.variable.updateValue (indigo.variables[76116845], str(Fcd)) # "cloudFactor"
   indigo.server.log("Updated Fcd = " + str(Fcd))
   if Fcd <= 0.05:
      Fcd = 0.05
      indigo.server.log("using lower limit value for  Fcd = " + str(Fcd))
   elif Fcd > 1: #This shouldn't happen, warn user if it does
      Fcd = 1 #actual radiation can't be greater than clear sky radiation
      indigo.server.log("WARNING - using upper limit value for  Fcd = " + str(Fcd))
   indigo.variable.updateValue (indigo.variables[76116845], str(Fcd)) # "cloudFactor"
else: # sun too low - use stored value for cloudiness factor
   indigo.server.log("using stored value for  Fcd = " + str(Fcd))   


# Calculate net long-wave radiation radiated upward using eq. 44
Rnl = sigma * Fcd * (0.34 - 0.14 * math.sqrt(ea)) * (avgTempK ** 4)
indigo.server.log("Rnl = " + str(Rnl) + " MJ/m2h-1")
Rns = Rs * (1- albedo) # correct shortwave radiation for albedo of reference crop
indigo.server.log("Rns = " + str(Rns) + " MJ/m2h-1")
Rn = Rns - Rnl # Calculate net radiation
indigo.server.log("Rn = " + str(Rn) + " MJ/m2h-1")

# Calculate Soil Heat Flux Density (G) using eqs. 65a, 65b (grass reference crop)
# Use different values depending on whether there was any solar radiation in the period
if Rn > 0:
   G = 0.1 * Rn
else:
   G = 0.5 * Rn

#And FINALLY we can calculate the evapotranspiration for this hour (ET) using eq. 1
if Rn > 0: # Use daytime value for Cd
   ET = (0.408 * delta * (Rn - G) + gamma * (Cn/(avgTempC + 273) * U2 * (es - ea))) / (delta + (gamma * (1 + Cd_day * U2)))
else: # Use nighttime value for Cd
   ET = (0.408 * delta * (Rn - G) + gamma * (Cn/(avgTempC + 273) * U2 * (es - ea))) / (delta + (gamma * (1 + Cd_night * U2)))
ET_in = (ET / 25.4)
indigo.server.log("Current hourly ETo (in mm) = " + str(ET))
indigo.server.log("Current hourly ETo (in inches) = " + str(ET_in))


# Update indigo variables
dailyETo = dailyETo + ET_in
runningETo = runningETo + ET_in
indigo.server.log("Current daily ETo (in inches) = " + str(dailyETo))
indigo.server.log("Current running total ETo (in inches) = " + str(runningETo))
indigo.variable.updateValue (indigo.variables[196835831], str(ET_in)) # "hourlyETo"
indigo.variable.updateValue (indigo.variables[1849325147], str(dailyETo)) # "dailyETo"
indigo.variable.updateValue (indigo.variables[1374037415],str(runningETo)) # "runningETo"

#  re-zero indigo running totals for observation period
indigo.variable.updateValue (indigo.variables[1514345125], str(0)) # "totalSolRad"
indigo.variable.updateValue (indigo.variables[900917121], str(0)) # "totalTempC"
indigo.variable.updateValue (indigo.variables[137383963], str(0)) # "totalHumid"
indigo.variable.updateValue (indigo.variables[1583976], str(0)) # "totalPress"
indigo.variable.updateValue (indigo.variables[27494067], str(0)) # "totalWind"
indigo.variable.updateValue (indigo.variables[582955114], str(0)) # "totalSampleCount"


You'll also need to totalize your hourly weather station data (or pull it from your weatherstation DB).

Comments, suggestions, bug reports welcomed.

...I opted to calculate the reference evapotranspiration constant based on hourly averages of weather data. I've had the thing working for a few days now, and it seems to agree with console calcs and common sense. If anyone who has a weather station that measures temperature, relative humidity, wind speed, and solar radiation wants to give it a try here it is...[snip]

WOW! I am totally impressed.

I can't tell you how many times I started to implement an ETo calculation on my own and just flailed. Finally, when I found the CIMIS data I figured why bother continuing to bang my head on the table and just used theirs. But, I'd really rather calculate my own -- the CIMIS data is roughly correct for my location, but is still a few km off. I do not have a weather station of my own, but I think it would be pretty easy to get the data required for your script from a nearby Weather Underground personal weather station (PWS).

Thanks for a GREAT contribution.

Thanks for this script!

I assume that all of the measurements coming from the Vantage are metric, not English or knots?

The secondary question -- I only have a Vantage Vue, so no access to solar radiation. I do know there are a couple of nearby Wunderground sites that are Vantage. I wonder if there's a way to use their solar radiation figure. Or maybe I could reverse calculate a solar radiation from my solar array's output in Watts... Hmmm.

You're welcome! Yes, I tried to document the units in the script, but everything is metric, it was easier that way.

Thanks!

I had a few minutes to look at PVOutput.org's API and I think I can do something with that data. I know next-to-nothing about scripting, so I'm sure this will be a painful exercise.

...[the] Davis Vantage Pro 2 weather station... ...console displayed evapotranspiration data, but WeatherSnoop (which is otherwise awesome) didn't support it... ...I opted to calculate the reference evapotranspiration constant based on hourly averages of weather data. I've had the thing working for a few days now, and it seems to agree with console calcs and common sense.

Just curious if your experience has changed since your first report. Do you continue to find a close tracking between the Davis ETo and HomeBrew calculations?

I was also wondering how you are adjusting your irrigation times based on the ETo. Specifically, what period ETo do you use and when do you apply it. For example, with the CIMIS data, early each morning I get their report for the previous 24 hours and then apply that against my irrigation times for that day. Essentially, rather than using the ETo factor for its own absolute value, I just trust that the relative changes on a day-to-day basis are correct. In the end, I do not know, or I guess really care, if the CIMIS data is for the preceding 24 hours, or just the daylight hours. However, with a local calculation updated hourly, I would have a lot more control over my choice of ETo period to use - and, in some cases, too much data may not be a good thing. Thus, I am very interested in hearing what you are doing and how it has ben working out.

You'll also need to totalize your hourly weather station data (or pull it from your weatherstation DB).

I assume you are using WeatherSnoop and the WeatherSnoop plugin to pull the data into Indigo variables. Can you confirm that.

Thanks.

You'll also need to totalize your hourly weather station data (or pull it from your weatherstation DB).

I assume you are using WeatherSnoop and the WeatherSnoop plugin to pull the data into Indigo variables. Can you confirm that.

Based on a hint in the script, that's what I wound up doing. The script pulls all of its data from Indigo variables, so if you have the ability to get the source data into Indigo, the script will work.

Basically, I used the WeatherSnoop plugin to get the current value and then used Python to add to existing value.

Code: Select all

totalTempC = int(float(indigo.variables[1925429375].value))
currentTempC = int(float(indigo.variables[1601247264].value))
cumulativeTempC = totalTempC + currentTempC
indigo.variable.updateValue(1925429375, value=str(cumulativeTempC))

This is going fairly smoothly. I have an ETo number being generated and with the WeatherSnoop plugin, I'm able to subtract out the rainfall since the last sprinkler cycle. I'd like to be able to subtract out the forecasted precipitation amount (not probability), but unless there's a way to get that from NOAA's XML feed, I'm not sure how I'd do that.

Otherwise, I'm in the process of adding zones and adjusting each zone for crop coefficient and microclimate. I had already calculated these from my WeatherTrak system.

Code: Select all

runningEToMM = int(float(indigo.variables[1098271889].value))
runningRain = int(float(indigo.variables[849525392].value))
totalIrrigationRequired = runningEToMM - runningRain
indigo.variable.updateValue(1437705013, value=str(totalIrrigationRequired))

z1kc = (float(indigo.variables[1583698197].value))
z1mc = (float(indigo.variables[491179291].value))
z1IrrigationRequired = totalIrrigationRequired * z1kc * z1mc
indigo.variable.updateValue(2132946, value=str(z1IrrigationRequired))

From there, the next step will be creating sprinkle and soak cycles for the zones and then I can tie it all together with some script to make only one zone work at a time.

Still need to pick an irrigation system that can handle 12 zones, work independently if necessary and also work with Indigo.

...Otherwise, I'm in the process of adding zones and adjusting each zone for crop coefficient and microclimate....
...From there, the next step will be creating sprinkle and soak cycles for the zones and then I can tie it all together with some script to make only one zone work at a time.

One note re: ETo control. I have not seen any systems that offer ETo control on a per-zone basis. In fact, as I noted in another post on the newAQUA system, it may be the only one that might, in the near future, offer it as a global setting.

What this (may) mean for your setup is that rather that getting a "smart" system land overriding it, you may be better off with a "dumber" system (like the EZFlora) and controlling it completely from Indigo or, possibly, the IrrigationCaddy - IF the Indigo plugin can be updated to support loading a complete schedule. For me, if the netAQUA can add a global duration multiplier, I will probably go that route and just hope that down the line they expand that to a per-zone multiplier. I know they are interested in getting WaterSense certification and I think that would require a fairly fine granularity of control.

The early adopter's mantra: One Can Hope

...Otherwise, I'm in the process of adding zones and adjusting each zone for crop coefficient and microclimate....
...From there, the next step will be creating sprinkle and soak cycles for the zones and then I can tie it all together with some script to make only one zone work at a time.

One note re: ETo control. I have not seen any systems that offer ETo control on a per-zone basis. In fact, as I noted in another post on the newAQUA system, it may be the only one that might, in the near future, offer it as a global setting.

The only system I know of that does this is from the WeatherTrak family (http://www.hydropoint.com) They download ETo nightly by radio and dynamically change the schedule for each of the zones based on the parameters you've configured. It's a great system, but it's very expensive and the residential versions are black boxes -- you really have zero control on what's going on -- definitely none if you're not in front of the box. And also no reporting, which is annoying.

The early adopter's mantra: One Can Hope

Indeed!

I will look again at EZFlora, but I would need two of them, which makes it more expensive, and it seems less flexible than some of the other options. At least I know Insteon well, so there would be few surprises in store.

I've completed a beta watering schedule in Indigo that takes the ETo, sprinkler mm/min, kc and microclimate information and applies it specifically to each zone. It's not as sophisticated as WeatherTrak, which also takes root depth and slope into their models, but at least I can estimate how much water I would use with my own program.

...For me, if the netAQUA can add a global duration multiplier, I will probably go that route and just hope that down the line they expand that to a per-zone multiplier. I know they are interested in getting WaterSense certification and I think that would require a fairly fine granularity of control...

I have been in contact with NetAQUA and can report they are planning to add support for additional weather-based control sometime this spring. It is not clear how much arbitrary external control will be available (I.e. some level of API). Although, if they do the weather right, the need for Indigo intervention is minimized, and the few necessary operations, like turn a zone program or schedule off can be accomplished through curl commands to the web interface.

After checking out all of the other recent posts on this thread, it still looks to me like NetAQUA has the most robust, well thought out, and finest grained irrigation control in the field -- If they can deliver the promised expansion capabilities and improved dynamic weather based controls I think it is a winner.

More news as it becomes available...

(CAVEAT: I am not looking for an irrigation system with a full-API. Rather, I am looking for a robust system which allows dynamic weather based irrigation that can be overridden and monitored by Indigo and my comments are made in that vein.)

Cycling back around to this topic as I make adjustments for drought and summer. I was wondering if anyone has compared the WeatherSnoop3 ETo data to the output from the script posted by @Grognard?