Calculate evapotranspiration

Question

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

Answer 1

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

Answer 2

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.

Answer 3

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)