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
The files in the water directory contain the amount of irrigation performed in mm. If the value is -1, the calculation is reset at that time. This is useful for start of the season when you don't care about past history anymore.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Thank you Dave. I look forward to checking those out. Your method to determine Solar Radiation is perfect.

I am just about to the point of comparing the figures to others in my area. Current / hourly Solar Radiation figures from within a few miles of here are available online from a nearby college.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Its working well Dave. I did look through your resources. You did an amazing job finding and interpreting those sources.

One question... If it rained say... 100 mm in a day. You wouldn't have to water the lawn again for a long time. Is that the theory behind ET? Or is that not practical. Is it just a matter of your soil's ability to absorb and hold the water?

I was thinking about decreasing the ending balance of older days by some sort of increasing percentage by their age. But that seems to circumvent the whole ET concept. What did you come up with in that regard?
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Good question! That is idea of soil saturation . At some point, the water just rubs off and cannot be absorbed anymore. We should add something to effectively 'cap' high amounts of rainfall. Up here in Ottawa, we should also try and manage spring thaw and the fact that the ground cannot absorb water since it is frozen in the spring. Both of these are handled manually wiry the '-1' reset.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
As I thought of it more last night, I too came up with the cap idea. I was thinking what if you had a statement that said something like if the rainfall was greater than three days of watering, then it would be equal to that equivalent. For example, if > 18 mm, then = 18. (6mm per day average????) And I don't know what to think about the frozen issue.

A person on the Homeseer board is interested too. I provided a link to this thread. He writes as though he has gone through the same issues we did. But I was lucky enough to find your thread. Thanks! I also told him I would post my adaptation of your script that includes links to Homeseer.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Yeah, maybe we could include a cap in the code that just says that if it ever rains more than X mm in one day (or per hour, in theory), that we cap the rainfall to be X.

The real value of X is dependent on the soil type in the area, and how fast the rail actually falls. More info here: http://www.fao.org/docrep/r4082e/r408...

For the frozen issue, I think I would just rely on using the reset. I'm not going to turn the water on until the ground is thawed anyway.

You can also try playing around with the "window" parameter to change how far back in time you look, but I don't think I would start decreasing the rainfall or solar radiation for days in the past, since it will give unpredictable results.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
I put a cap of 18 mm per day. The balances look a little better that way. I messed with the window value quite a bit but if it rained 45 mm yesterday, the window doesn't help.

I have it set up such that you can select any window you please and the script will set the previous day's values to zeros. Maybe the window should just be kept at a constant of ...... 5 days? 4 days? And never look further back than that? And the values previous to the start date would be zeros?

Also someone asked about elevation. I see from your documentation - a table in annex 2 - that adjusts for various altitudes other than the default of 100 m. I am wondering whether that parameter should be something you input.

I hope the grass appreciates this effort.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
From what I could tell, the main effect of elevation is the pressure change, and I think I take that into account with the "g = 0.665e-3 * pressure", but somebody can check my math.

There are assumption and simplifications in the implementation. If you have the actual data, you can look in www.fao.org/docrep/x0490e/x0490e07.htm to try and be more accurate.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Hi Dave,

I compared 40 days of your readings with those of a college weather station. The comparison I used involves two locations about 8 miles apart. I can improve that by choosing a weather station that is closer to the college station (college station posts hourly updates of their readings of ET and SR).
http://mesonet.k-state.edu/weather/hi... . This is a station in Olathe, Kansas. I should have used a weather station in Olathe as well, which I will do.

The Solar Radiation results are amazingly close. I was surprised. The ET figures aren't quite as close. I need to compare your readings with an online ET calculator maybe. Have you ever done that? Perhaps I've made a mistake or perhaps something is not quite right with the ET formula. I don't know what could be different about the 'type' of ET reading, but I don't know enough about it to have an idea. I multiplied yours by .039370 to get inches.

I will download 60 days of data to compare from a wx closer to the college station.

Getting the ET should be easy as compared to getting the SR, which you have done very well.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Thanks! I will take a closer look. I could very well have a bug. It seemed that the et numbers were a little high to me, but didn't have any way to compare.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
Sounds good. As a very rough estimate, it seems if your ET figures were multiplied by .7, the average difference over a 158 day comparison would be 0%.

I downloaded data from the K-State wx station back to May 1.

But I need to more closely review the differences. And I thought I would figure out a way to substitute their SR figures for yours and see if your ET figures more closely track theirs even if they are high. If that worked out and you found a bug with the ET calculation, then we could go back and look to tweaking the SR. But overall, your SR looks good.

Happy Hunting
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
I've been through everything in your formulas except for the solar radiation section and can't find anything wrong.

There are some fixed parameters that I'm not sure how to deal with. One is the equation 38 in table 2.8. Do you know how to implement that equation?

Did you consider doing this whole thing using the hourly process rather than daily? You have the foundation built for hourly with your cloud cover calculations. The reason I ask is that the K-State solar radiation and ET data is available on an hourly basis. Perhaps using that would point to the actual difference between theirs and yours.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
How far off is it? Can you post some comparisons so I can get a better feel for what might be missing?

It looked like the right ballpark to me.
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Equation 38 looks like it is accounting for roughly 23% of the incoming solar radiation to be reflected by the grass, which in theory could decrease the net solar radiation by be 77%. There is also a net longwave radiation that seems to be emitted. You can try emailing me some data directly. My gmail account name is the same as my wunderground name

It looks like I am using the Rn wrong. I calculate the solar radiation taking into account cloud cover, but Rn is the net longwave radiation. It looks like they may ignore that in some cases, but it could be the main difference you are looking for.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
That would be great if you found something wrong.

I just emailed three Excel spreadsheets to you. If you don't have Excel, let me know.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
emailed one more sheet that shows the result of your ET readings reduced to 70% of their values
Photo of davemiedema

davemiedema

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

This gist contains updated equations based on some digging from Frank (thanks!). It looks like better results than we were getting before.
Photo of David Achterberg

David Achterberg

  • 3 Posts
  • 0 Reply Likes
I found your script yesterday while searching for something similar and downloaded the version at the github site. I am running Yosemite on a Mac and testing with the terminal using Python 2.7.6. I was going through it last night and corrected a few formatting issues in bringing it over, but I have a few comments/questions. I am not proficient in Python, so I hope these are not too obvious.

I assumed the references to the api needed to be changed to the imperial versions for the US, which I did.

I changed the api url in the script to forecastURL = 'http://api.wunderground.com/api/' + KEY + '/forecast10day/q/' + str(ZIPCODE) + '.json', since the version in the script didn't seem to get any location data. I left the history url (further along in the script) as is, it seems to get the history when used in a browser.

I commented out the following lines, since I need to investigate this more. I know from reading the post that the script will pull from and write to a local file, but I don't understand the formatting and type of the file. Do you have an example? Anyway, commenting it out for now doesn't error out, but I'm not sure if this is the correct way to bypass it for now.
# response = urllib.urlopen(historyURL).read()
# cachefile = open("data/" + datestring, 'w')
# cachefile.write(response)
# cachefile.close()

I get the error:
in getHistoricalData    if (cloudCover == -1): cloudCover = previousCloudCover
UnboundLocalError: local variable 'previousCloudCover' referenced before assignment
I changed the line to read previousCloudCover = 0, which is not right, but it seems that it could loop with previousCloudCover = -1, and changing to 0 doesn't error out.

I get the error:
in getHistoricalData    Rnl = aveSigmaT * (0.34 - 0.14 * math.sqrt(ea)) * (1.35 * Rs / Rso - 0.35)
ZeroDivisionError: float division by zero
This I have no clue on how to fix, and was hoping for some direction. It seems necessary in the calculation of ET and without going to the websites you referenced I was hoping for a simple fix. I think if I get over this, I can get the script to complete.

Sorry for the verbosity, and that was a few more than I intended to bring up, but I wanted to let you know what issues I was getting. Since I was calling the api anyway, I am trying to incorporate another script in the def getForecastData() function that will update a variable based on the rain forecast/history that is used in another calculation to determine runtime for an irrigation system.

Thanks for your work on this as well.
Photo of frankpc

frankpc

  • 25 Posts
  • 0 Reply Likes
I learned Python just to work on this.  I have since moved on to visual basic, which is more compatible with Homeseer, which is what I use to control the script.  Before I respond to this, do you have an interest in the visual basic version?  It is much more complete and provides more accurate results.  In addition, it figures ET, RS, and RSo, etc., several different ways.  I believe it also has more comments.

I will email the script to you if I had your email address.

The original written by Dave is what got several people involved with ET.  I too thank him for his work on this.  His script is what inspired the rest of us to get involved with ET.
(Edited)
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
All the units need to be metric.

Commenting out this

     response = urllib.urlopen(historyURL).read()

will mean that the history API will never be called, so you won't get any historical data.  That likely explains the divide by zero you are getting.  It probably also explains the use of an uninitialized previousCloudCover, since there is no data to fall back to.

You will need to fix the URLs to remove references to Ottawa, and put your own URL references in for wherever you want the data from.  I'm not sure why I didn't use the COUNTRY and CITY variables in the historyURL, but I should have.

Also, I am not sure if the data returned in forecast10day API is the same as the forecast API.
Photo of David Achterberg

David Achterberg

  • 3 Posts
  • 0 Reply Likes
Thanks for the response.
I changed all the units back to metric (makes sense now reading through the calculations); even changed ELEVATION to meters. I don't know if GMTOFFSET is positive or negative? It should be -6 for me, but i left it positive. I tried running the script both ways.
I didn't note before but I left the line "data = json.loads(urllib.urlopen(historyURL).read())" uncommented, that should get the data anyway, correct? But irrelevant now, since I fixed the "cachfile = open syntax" (another comment below).
I changed "forecastURL = 'http://api.wunderground.com/api/' + KEY + '/forecast/q/' + str(ZIPCODE) + '.json'" instead of the 10 day forecast and it pulls in a slightly different format, but works, since "forecastURL = 'http://api.wunderground.com/api/' + KEY + '/forecast/q/' + COUNTRY + '/' + CITY + '.json'" doesn't pull in any forecast for me.
I restored the original previousCloudCover syntax.
I changed the "cachefile = open" syntax to create/write to a text file (I had the path syntax wonky before), and the daily history gets written to it now.

So all of the original script should be restored. The error I get now is "in getHistoricalData
    solarRadiation = clearSkyInsolation * (1 - 0.75*(math.pow(cloudCover,3.4)))
ValueError: math domain error"

I'm still doing something wrong surely, but I double checked the code. Would you have any ideas on where to look?

Thanks again for the help.
Photo of David Achterberg

David Achterberg

  • 3 Posts
  • 0 Reply Likes
frankpc:
Thanks for the offer, but I'm not sure visual basic will work on a mac? Although I shouldn't look a gift horse in the mouth, so to speak. I'm trying to get up to speed with python, so I would like to stick to it for the moment. I wouldn't mind looking at what you have anyway, but does this forum have PM to get you an email address?
Photo of davemiedema

davemiedema

  • 29 Posts
  • 7 Reply Likes
Not sure.  Can you print out the values of all the parameters to that equation and see if one looks wrong?  math domain errors are usually when you do an undefined math operation (square root of a negative number, or log of a negative number).
Photo of wingdspur

wingdspur

  • 4 Posts
  • 1 Reply Like
Here are two libraries, one for Python and the other is a Ruby port of it for calculating reference crop evapotranspiration (ETo), sometimes referred to as potential evapotranspiration (PET). The library provides numerous functions for estimating missing meteorological data.

Three methods for estimating ETo/PET are implemented:

  • FAO-56 Penman-Monteith (Allen et al, 1998)
  • Hargreaves (Hargreaves and Samani, 1982; 1985)
  • Thornthwaite (Thornthwaite, 1948)
You can estimate solar radiation values with functions in these libraries.

Python - PyETo: https://github.com/woodcrafty/PyETo

Ruby - evapotranspiration: https://github.com/brycejohnston/evapotranspiration
Documentation http://www.rubydoc.info/gems/evapotranspiration/Evapotranspiration/FAO
(Edited)