Calculating degree days with NumPy

Michael Roberts offers an extremely helpful post at G-FEED, including R code for generating the season growing degree days measure used in the seminal Schlenker Roberts 2009 paper on nonlinear crop responses to temperature. Since many of the papers claiming to follow the approach don’t actually adhere to it, he’s posted code in the interest of encouraging people to use the sinusoidal fitting method and perhaps improve on their standard errors.

I’ve been reliant on the wonderful Climate Data Operators library for working with netCDF files and have therefore been basing most of my current weather/climate data analysis in Python (see the Python bindings for CDO here). I thought it might be helpful to provide the Python analogue of the R code so that folks in a Python workflow might be up and running. I’ve modified the function to address a single pixel at a time (i.e., you’ll need to write the loop that spans all relevant grid cells), and hence there’s no cell count variable here to average over.

Head over to Michael’s original post, since I’ve offered the code snippet below without documentation. Drop a comment if you find this to be a helpful addition.

[code language=”python”]
def days_in_range(t0, t1, tMin, tMax):

n = len(tMin)
t0 = np.array([t0] * n)
t1 = np.array([t1] * n)
t0[np.where(t0 < tMin)] = tMin[np.where(t0 < tMin)]
t1[np.where(t1 > tMax)] = tMax[np.where(t1 > tMax)]
outside = (t0 > tMax) | (t1 < tMin)
inside = np.logical_not(outside)

def u(z, ind, t_min, t_max):
z_t = np.array([val for x, val in zip(ind, z) if x])
tMax_t = np.array([val for x, val in zip(ind, t_max) if x])
tMin_t = np.array([val for x, val in zip(ind, t_min) if x])
return np.divide(z_t – tMin_t, tMax_t – tMin_t)

time_at_range = (2/pi)*((np.arcsin(u(t1, inside, tMin, tMax)) – np.arcsin(u(t0, inside, tMin, tMax))))
return np.sum(time_at_range)

[/code]

To ensure that all results from calculations in the functions (e.g., the vector division in the u function) produces floats, ensure that

[code language=”python”] t0, t1 [/code]

are treated as floats by including a decimal afterwards (the range of 21-22C should be 21.0 and 22.0).

In a follow-up post, I’ll include snippets that span all pixels in the dataset and summarize the data and folder management system I’ve been using when exporting GDDs for use by Stata or ArcGIS.

I disclaim any responsibility for issues (emotional, financial, or other) associated with use of this code.