Calculate evapotranspiration

I am trying to create a code for the Basic evapotranspiration equation by Hargreaves. I have attached a screenshot of the equation to be replicated.

I want to create a class wherein I can input data which consists of tmin, tmax, tmean, day and lat. I have written the code below based on what I can understand from the equation.

import math
import numpy as np

class Hargreaves:

#Instantiating the class

def __init__(self, tmin, tmax, tmean, day, lat):
self.tmin = tmin
self.tmax = tmax
self.tmean = tmean
self.day = day
self.lat = lat

# latrad = pyeto.deg2rad(lat)  # Convert latitude to radians
# day_of_year = datetime.date(2014, 2, 1).timetuple().tm_yday

# setting the gsc value

@property
def gsc(self):
return 0.082

# Calculating d

@property
def d(self):
return 0.409 * math.sin((2*math.pi) * (self.day/365) - 1.39)

# Calculating ws

@property
def ws(self):
return math.acos(-math.tan(self.lat) * math.tan(self.d))

# Calculating dr

@property
def dr(self):
return 1 + 0.033 * math.cos((2*math.pi) * (self.day/365))

# Calculating Radiation

@property
def ra(self):
return 24 * 60 * (self.gsc/math.pi) * self.dr * self.ws * math.sin(self.lat) * math.sin(self.d) + (math.cos(self.lat) * math.cos(self.d) * math.sin(self.ws))

# Function to calculate evapotranspiration

@property
def et(self):
return (0.0023/2.45) * (self.tmean + 17.8) * ((self.tmax - self.tmin) ** 0.5 ) * self.ra

• I have concerns about the accuracy of that math (fixed year durations?), but that's a separable issue and I'll have to assume that you need an implementation faithful to this reference. Oct 5, 2021 at 3:55
• @Reinderien Correct - this existing implementation the author of the question seems to be aware of (some code copied verbatim!) does have proper conversion to Julian Days (the now-commented .tm_yday)... Oct 5, 2021 at 8:31
• @Reinderien Yes, I had previously converted the day to Julian days, but after looking at my data source, I found that the day variable will be provided in Julian days. Oct 5, 2021 at 9:06
• Why does math in your image seem so complicated but in code so easy? xD Oct 5, 2021 at 9:50
• Units are wrong for solar constant and extraterrestrial radiation. a/b/c is not the same as a/(b/c). Oct 6, 2021 at 19:36

3 Answers

I wonder if there is a need for a class here? It seems class Hargreaves has no responsibilities and is solely used to compute the evapotranspiration value based on the inputs provided in the constructor. Should a function be used instead?

import math
import numpy as np

def hargreaves(tmin, tmax, tmean, day, lat):
gsc = 0.082
dr  = 1 + 0.033 * math.cos((2*math.pi) * (day/365))
d   = 0.409 * math.sin((2*math.pi) * (day/365) - 1.39)
ws  = math.acos(-math.tan(lat) * math.tan(d))
ra  = 24 * 60 * (gsc/math.pi) * dr * ws * math.sin(lat) * math.sin(d) + (math.cos(lat) * math.cos(d) * math.sin(ws))
et  = (0.0023/2.45) * (tmean + 17.8) * ((tmax - tmin) ** 0.5 ) * ra
return et

• I'd agree if it were the only equation being so treated. In a complex system I'd however expect an interface "Equation" with a calculate method that takes another interface "Parameters" as input and returns an interface "Output" or something of the kind, then multiple implementations of that that can be plugged in as needed. But maybe I'm overthinking things. Oct 5, 2021 at 6:11
• Great solution, but I want to try to implementing a class for this particular equation. Oct 5, 2021 at 9:07
• @jwenting Given that functions are first class types in python I suspect you are overengineering. It sounds like you are trying to duplicate <class 'function'>
– Pace
Oct 6, 2021 at 5:04

To add to Peter's excellent answer, please also consider the likely audience of your code. I do not know who Hargreaves is, (or even that it's a name, looking only at the source), and practically every single name in the overall equation is completely meaningless to me.

If the only people who will ever read your code are going to be other people with domain knowledge sufficient to understand what, for example, gsc is at a glance, then no change is needed. If someone without domain knowledge might ever read it, they're going to see a bunch of impenetrable nonsense filled with "bad" variable names and magic numbers.

In the latter case, consider adding a docstring (preferred) or a comment with a link to some resource from which a completely uninitiated reader can understand what they need to. Such a reader might be a team member at work, assigned to check your code for bugs - they'd need to understand what proper operation of the code is meant to look like.

• This is typical with science code; I've written a lot of functions like that when I was a researcher. Sometimes descriptive names are better, but like you said, sometimes it's better to add a reference like "Variable names as in [paper citation] page N, equations X to Y", since anyone who needs to mess with the implementation will almost certainly also look up the paper the code is based on. In that case the code can be easier to read if it matches the spec, rather than having to translate concepts. Oct 5, 2021 at 7:49
• I'm having trouble finding something actionable in your comment. It seems to only affirm what I posted. Oct 6, 2021 at 8:58
• It was aimed at readers of your answer, with the intention to underline and expand a little on the rationale behind doing what you suggested. Oct 6, 2021 at 9:16
• @Dronir although it's "typical with science code" i do believe it's not at all recommended. it "matches with the spec" however unless the spec gives pseudocode which the code is a direct port of, generally i believe it just reads better with full names (main upside being not having to memorize the specific abbreviations used). it's also just good convention in general to name variables properly; naming them badly leads to bad habits, potentially leading to passing said bad convention on when teaching. Oct 6, 2021 at 22:42
• very important also to note that for this question specifically the full names are right next to the abbreviations. so the only potential downside is that the lines become unreadable because they are too long. Oct 6, 2021 at 22:42

As someone who has to read/write science code, please please include units in your constants, and in the output from your function. They are in the equations, but not the code.

I still wish that we had (in major circulation) languages that enforced unit checks! Would have avoided an expensive space crash, at the very least (Mars climate orbiter)

• Better yet, use a measurement abstraction which ensures correctness. e.g. python-measurement.readthedocs.io/en/latest/topics/use.html Oct 5, 2021 at 16:16
• python-measurement is not really ideal. w=Weight(kg=1) then w^0.5 yields an error. Basic operators don't work, and are not checked. Similarly math.ceil, fabs and friends don't work on such objects. Great where they work, otherwise just cumbersome. Oct 6, 2021 at 11:59
• That was just a quick example I googled, I don't normally work in Python. The language does support operator overloading, so presumably you can make a library that lets you treat these objects like numbers. Oct 6, 2021 at 12:03
• I would suggest Pint pint.readthedocs.io/en/stable to deal with units Oct 6, 2021 at 13:04
• "please please include units in your constants". Yes, especially since the units are wrong in the specifications: MJ/(m²/min) should be MJ/(m² * min). Oct 6, 2021 at 20:40