8
\$\begingroup\$

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. Hargreaves equation

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
\$\endgroup\$
6
  • 3
    \$\begingroup\$ 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. \$\endgroup\$
    – Reinderien
    Oct 5 at 3:55
  • 1
    \$\begingroup\$ @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)... \$\endgroup\$
    – Lodinn
    Oct 5 at 8:31
  • \$\begingroup\$ @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. \$\endgroup\$ Oct 5 at 9:06
  • \$\begingroup\$ Why does math in your image seem so complicated but in code so easy? xD \$\endgroup\$
    – Allart
    Oct 5 at 9:50
  • 1
    \$\begingroup\$ Units are wrong for solar constant and extraterrestrial radiation. a/b/c is not the same as a/(b/c). \$\endgroup\$ Oct 6 at 19:36
35
\$\begingroup\$

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
\$\endgroup\$
3
  • 2
    \$\begingroup\$ 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. \$\endgroup\$
    – jwenting
    Oct 5 at 6:11
  • \$\begingroup\$ Great solution, but I want to try to implementing a class for this particular equation. \$\endgroup\$ Oct 5 at 9:07
  • 2
    \$\begingroup\$ @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'> \$\endgroup\$
    – Pace
    Oct 6 at 5:04
10
\$\begingroup\$

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.

\$\endgroup\$
5
  • 10
    \$\begingroup\$ 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. \$\endgroup\$
    – Dronir
    Oct 5 at 7:49
  • \$\begingroup\$ I'm having trouble finding something actionable in your comment. It seems to only affirm what I posted. \$\endgroup\$ Oct 6 at 8:58
  • \$\begingroup\$ 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. \$\endgroup\$
    – Dronir
    Oct 6 at 9:16
  • \$\begingroup\$ @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. \$\endgroup\$
    – somebody
    Oct 6 at 22:42
  • \$\begingroup\$ 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. \$\endgroup\$
    – somebody
    Oct 6 at 22:42
7
\$\begingroup\$

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)

\$\endgroup\$
6
  • 2
    \$\begingroup\$ Better yet, use a measurement abstraction which ensures correctness. e.g. python-measurement.readthedocs.io/en/latest/topics/use.html \$\endgroup\$
    – Alexander
    Oct 5 at 16:16
  • \$\begingroup\$ 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. \$\endgroup\$ Oct 6 at 11:59
  • \$\begingroup\$ 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. \$\endgroup\$
    – Alexander
    Oct 6 at 12:03
  • 1
    \$\begingroup\$ I would suggest Pint pint.readthedocs.io/en/stable to deal with units \$\endgroup\$ Oct 6 at 13:04
  • 1
    \$\begingroup\$ "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). \$\endgroup\$ Oct 6 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.