Calculating degree days with NumPy

2010-08-08 at 11-41-43

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.

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)    

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

 t0, t1 

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.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s