6.6
"dd_table" Script
#!/usr/bin/python
# * coding: utf8 *
# dd_table.py
# 2018 Feb 23 . ccr
u"""Print Growing DegreeDay Table.
This script generates a table of growing degreedays used for modeling
the emergence of larvae of Cydia pomonella (codling moth).
"""
# 2019 Jan 08 . ccr . Change function names. Provide
# . . gdd_single_sine_vertical_cutoff entry point
# . . among others.
# 2018 May 28 . ccr . Compare degreedays of cooling.
from __future__ import division
import sys
import codecs
import math
ZERO = 0
SPACE = ' '
NULL = ''
NUL = '\x00'
NA = 1
STDOUT = codecs.getwriter('utf8')(sys.stdout)
FMT = '%3i'
def gdd_single_sine_with_theta( # 2019 Jan 08
day_max_temp,
day_min_temp,
threshold_temp,
):
u"""Return growing degreedays (GDD) above threshold, given daily max
and min temps.
This function is transcribed from a *perl* subroutine:
o Coop, Len. "DegreeDay Calculation Formulas." 8 May
2015. Integrated Plant Protection Center, Oregon State U. 27
Feb. 2018 <http://uspest.org/wea/wea3.html>.
This function has been verified against various published tables of
developmental degreedays for Cydia pomonella. It assumes that
diurnal temperatures follow a single symmetric sine wave. The
document (above) references this seminal article:
o Baskerville, G.L., and P. Emin. "Rapid Estimation of Heat
Accumulation from Maximum and Minimum Temperatures." _Ecology_
50 (1969): 514517.
I wouldn't say 1969 was early days for generalpurpose computing,
but apparently the journal editors for _Ecology_ were way more
interested in graphs and numerical tables than in computation
algorithms. Although the authors (B&E) mention a FORTRAN program,
it is not reproduced in the body of the article.
Here's how it works, though: The sine wave is clipped horizontally
from below by the threshold temperature. Growing degreedays for
each day is proportional to the area under the curve above the
threshold, which thus depends on the threshold and the maximum and
minimum temperatures for each day.
The algebra, trigonometry, and Integral Calculus are apparently not
too exotic but are fairly obscure. Here follows my reconstruction.
=====
If the threshold temp is BELOW the daily min temp, the answer we
must give is the conventional one:
dd_conventional = average_temp  threshold
where average_temp = (max + min) / 2
... because the curve is not clipped.
Below, we're going to use the following ¨simplifying¨ expression, so
I suppose I ought to define it here.
Let d2 = 2 * dd_conventional =
2 * threshold  (max + min)
Consider a sine wave from π/2 to 3π/2 radians with its peak at high
noon (π/2 radians).
If the threshold temp is ABOVE the daily min temp, we need to
integrate over a portion of the curve. The essential magic is to
discover the angle θ where the sine wave, scaled and offset to the
average temperature, crosses this threshold. Then we can cast away
the area below the threshold along with the tails of the curve.
Because the sine wave is symmetric about π/2, we need to integrate
only between π/2 and π/2.
The definite integral of the sine wave from π/2 to +π/2 is ZERO
because half of it is negative, so that doesn't do us much good.
We'll have to choose a different function. How about sine + 1?
The definite integral of (sine + 1) from π/2 to +π/2 is:
cos(π/2)  (cos(π/2)) + π/2  (π/2) = π
This is the area under the curve. To normalize the calculations to
follow, we have to multiply by this:
downscale_factor = 1 / π
When threshold temp == min temp and θ == π/2, the answer we must
give is dd_conventional. This becomes our upscale_factor:
upscale_factor = ((max + min) / 2)  threshold =
(max + min  2 * threshold) / 2 =
(max + min  2 * min) / 2 =
(max  min) / 2
The following formula is similar to that given by B&E in their
Fig. 1.B:
heat_units = downscale_factor * (upscale_factor * whole_area  clipped_area))
where:
whole_area = integral(dt) of (sin(t) + 1) from θ to π/2 =
(cos(π/2)  (cos(θ)) + (π/2  θ)) =
(cos(θ) + (π/2  θ))
clipped_area = integral(dt) of (threshold  min) from θ to π/2 =
(threshold  min) * (π/2  θ)
thus:
heat_units =
downscale_factor * (upscale_factor * (cos(θ) + (π/2  θ)))  (threshold  min) * (π/2  θ) =
(1 / π) * (((max  min) * (cos(θ) + (π/2  θ)) / 2)  ((threshold  min) * (π/2  θ))) =
(1 / 2π) * ((max  min) * (cos(θ) + (π/2  θ))  2 * (threshold  min) * (π/2  θ)) =
(1 / 2π) * ((max  min) * cos(θ) + (max  min) * (π/2  θ)  2 * (threshold  min) * (π/2  θ)) =
(1 / 2π) * ((max  min) * cos(θ) + ((max  min)  2 * (threshold  min)) * (π/2  θ)) =
(1 / 2π) * ((max  min) * cos(θ) + (max  min  2 * threshold + 2 * min) * (π/2  θ)) =
(1 / 2π) * ((max  min) * cos(θ) + (max  2 * threshold + min) * (π/2  θ)) =
(1 / 2π) * ((max  min) * cos(θ)  (2 * threshold  (max + min)) * (π/2  θ)) =
(1 / 2π) * ((max  min) * cos(θ)  d2 * (π/2  θ))
... which is the result calculated by this function!
=====
QED
=====
Now, for the essential magic we need to compute:
θ = arcsin(unit_base)
... where unit_base is the threshold temp scaled to the sine wave
(not sine + 1):
unit_base = 2 * (threshold  min) / (max  min)  1 =
(2 * (threshold  min)  (max  min)) / (max  min) =
(2 * threshold  2 * min  max + min) / (max  min) =
(2 * threshold  min  max) / (max  min) =
(2 * threshold  (max + min)) / (max  min) =
d2 / (max  min)
But there's a fly in this ointment. We want to have pretty good
computational accuracy where unit_base is pretty close to ±1, which
should yield an arcsin result pretty close to ±π/2, but that is just
where the accuracy of realworld arcsin implementations falls down
due to inherent granularity of simulating mathematics of continuous
functions on binary machines, which can handle only discrete values.
Maybe it would be better to calculate θ using the arctan, instead.
I think it's just this easy:
sin(θ) = a / c where the sides of the right triangle are:
(a * a) + (b * b) = (c * c)
and, using the same angle:
tan(θ) = a / b = a / sqrt ((c * c)  (a * a))
Now θ = arcsin(a / c)
and, as above, this means:
a = unit_base and c = 1
Then, the same angle would be:
θ = arctan(a / sqrt((c * c)  (a * a)))
θ = arctan(unit_base / sqrt((1 * 1)  (unit_base * unit_base))) =
arctan((d2 / (max  min)) / sqrt(1  (d2 / (max  min)) * (d2 / (max  min)))) =
arctan(d2 / sqrt((max  min) * (max  min)  d2 * d2))
... which, if I have not deluded myself, is the formula for θ used
by the *perl* function referenced above and is thus the formula
that this script uses, too. (In practice however, using the arcsin
yields identical results rounded to the nearest integer.)
=====
Does it need to be this complicated?
It is unreasonable in cooler weather early in the season to require
that the AVERAGE temperature be above the threshold before recording
any degreeday accumulation because plainly the actual temperature
may actually be ABOVE the threshold for part of the day even when the
AVERAGE is less than the threshold. The single sine wave model
captures this small amount of heating whereas the conventional model
does not. Thus, in the earlier (and the later) portions of the
season the cumulative growing degreedays will (and should) outrun
the conventional degreedays.
Does it need to be any more complicated?
The actual GDD could be calculated, given hourly readings, but that
would be a lot more effort than just collecting the daily max and
min temps and considering each 24hour (midnight to midnight) day to
be independent of preceding and following days. This seems to be
accurate enough to predict insect life stages to within a few days.
An alternative would be to collect max/min readings for 12hour
periods.
"""
day_sum_temp = day_max_temp + day_min_temp
day_diff_temp = day_max_temp  day_min_temp
theta = None
if day_diff_temp < ZERO:
result = None
elif threshold_temp < day_min_temp:
result = (day_sum_temp / 2.0)  threshold_temp # dd_conventional
elif threshold_temp > day_max_temp:
result = ZERO
else:
d2 = threshold_temp + threshold_temp  day_sum_temp # 2.0 * dd_conventional
theta = math.atan2(d2, math.sqrt(day_diff_temp * day_diff_temp  d2 * d2))
if (d2 < ZERO ) and (theta > ZERO):
theta = theta  math.pi
# if day_diff_temp == ZERO:
# theta = ZERO
# else:
# theta = math.asin(d2 / day_diff_temp)
result = ((day_diff_temp * math.cos(theta))  (d2 * ((math.pi / 2.0)  theta))) / (math.pi * 2.0)
return (result, theta)
def gdd_single_sine_no_cutoff( # 2019 Jan 08
day_max_temp,
day_min_temp,
threshold_temp,
cutoff_temp=None,
):
(result, theta) = gdd_single_sine_with_theta( # 2019 Jan 08
day_max_temp=day_max_temp,
day_min_temp=day_min_temp,
threshold_temp=threshold_temp,
)
return result
def dd_conventional( # 2018 Jan 08
day_max_temp,
day_min_temp,
threshold_temp,
cutoff_temp=None,
):
u"""Calculate degreedays of cooling (DD).
This is a heating, ventilation, and airconditioning (HVAC) rule of
thumb and is often held out as an adequate and easily grasped
substitute for calculating growing degreedays.
"""
day_sum_temp = day_max_temp + day_min_temp
day_diff_temp = day_max_temp  day_min_temp
if day_diff_temp < ZERO:
result = None
else:
result = (day_sum_temp / 2.0)  threshold_temp
if result < ZERO:
result = ZERO
return result
def gdd_single_sine_horizontal_cutoff( # 2019 Jan 08
day_max_temp,
day_min_temp,
threshold_temp,
cutoff_temp,
):
u"""Implement both lower and upper threshold temps.
Test the gdd_single_sine_horizontal_cutoff function.
>>> int(gdd_single_sine_horizontal_cutoff(50, 48, 50, 88 ) + 0.5)
0
>>> int(gdd_single_sine_horizontal_cutoff(50, 50, 50, 88 ) + 0.5)
0
>>> gdd_single_sine_horizontal_cutoff(50, 52, 50, 88) is None
True
>>> int(gdd_single_sine_horizontal_cutoff(52, 48, 50, 88 ) + 0.5)
1
>>> int(gdd_single_sine_horizontal_cutoff(52, 48, 50, 88 ) + 0.5)
1
>>> int(gdd_single_sine_horizontal_cutoff(55, 48, 50, 88 ) + 0.5)
2
>>> int(gdd_single_sine_horizontal_cutoff(52, 51, 50, 88 ) + 0.5)
2
>>> int(gdd_single_sine_horizontal_cutoff(55, 51, 50, 88 ) + 0.5)
3
>>> int(gdd_single_sine_horizontal_cutoff(52, 52, 50, 88 ) + 0.5)
2
>>> int(gdd_single_sine_horizontal_cutoff(54, 48, 50, 88 ) + 0.5)
2
>>> int(gdd_single_sine_horizontal_cutoff(54, 50, 50, 88 ) + 0.5)
2
>>> int(gdd_single_sine_horizontal_cutoff(54, 52, 50, 88 ) + 0.5)
3
>>> int(gdd_single_sine_horizontal_cutoff(86, 48, 50, 88 ) + 0.5)
17
>>> int(gdd_single_sine_horizontal_cutoff(86, 50, 50, 88 ) + 0.5)
18
>>> int(gdd_single_sine_horizontal_cutoff(86, 52, 50, 88 ) + 0.5)
19
>>> int(gdd_single_sine_horizontal_cutoff(88, 48, 50, 88 ) + 0.5)
18
>>> int(gdd_single_sine_horizontal_cutoff(91, 48, 50, 88 ) + 0.5)
19
>>> int(gdd_single_sine_horizontal_cutoff(91, 51, 50, 88 ) + 0.5)
21
>>> int(gdd_single_sine_horizontal_cutoff(88, 50, 50, 88 ) + 0.5)
19
>>> int(gdd_single_sine_horizontal_cutoff(88, 51, 50, 88 ) + 0.5)
20
>>> int(gdd_single_sine_horizontal_cutoff(88, 52, 50, 88 ) + 0.5)
20
>>> int(gdd_single_sine_horizontal_cutoff(90, 48, 50, 88 ) + 0.5)
19
>>> int(gdd_single_sine_horizontal_cutoff(90, 50, 50, 88 ) + 0.5)
20
>>> int(gdd_single_sine_horizontal_cutoff(90, 52, 50, 88 ) + 0.5)
21
>>> int(gdd_single_sine_horizontal_cutoff(86, 86, 50, 88 ) + 0.5)
36
>>> gdd_single_sine_horizontal_cutoff(86, 88, 50, 88) is None
True
>>> gdd_single_sine_horizontal_cutoff(86, 90, 50, 88) is None
True
>>> int(gdd_single_sine_horizontal_cutoff(88, 86, 50, 88 ) + 0.5)
37
>>> int(gdd_single_sine_horizontal_cutoff(88, 88, 50, 88 ) + 0.5)
38
>>> gdd_single_sine_horizontal_cutoff(88, 90, 50, 88) is None
True
>>> int(gdd_single_sine_horizontal_cutoff(90, 86, 50, 88 ) + 0.5)
37
>>> int(gdd_single_sine_horizontal_cutoff(90, 88, 50, 88 ) + 0.5)
38
>>> int(gdd_single_sine_horizontal_cutoff(90, 90, 50, 88 ) + 0.5)
38
"""
dd_threshold = gdd_single_sine_no_cutoff(day_max_temp, day_min_temp, threshold_temp)
dd_cutoff = gdd_single_sine_no_cutoff(day_max_temp, day_min_temp, cutoff_temp)
if None in [dd_threshold, dd_cutoff]:
result = None
else:
result = dd_threshold  dd_cutoff
if result < ZERO:
result = ZERO
return result
def gdd_single_sine_vertical_cutoff( # 2018 Jan 09
day_max_temp,
day_min_temp,
threshold_temp,
cutoff_temp,
):
u"""Implement a vertical cutoff.
"""
(dd_threshold, theta_1) = gdd_single_sine_with_theta(day_max_temp, day_min_temp, threshold_temp)
(dd_cutoff, theta_2) = gdd_single_sine_with_theta(day_max_temp, day_min_temp, cutoff_temp)
if None in [dd_threshold, dd_cutoff]:
result = None
else:
result = dd_threshold  dd_cutoff
if (day_max_temp < cutoff_temp):
pass
elif (cutoff_temp < day_min_temp):
result = ZERO
else:
# integral(dt) of (cutoff  min) from θ to π/2
rectangle = (math.pi / 2.0  theta_2) * (cutoff_temp  threshold_temp) / math.pi
result = rectangle
if result < ZERO:
result = ZERO
return result
def gdd_single_sine_intermediate_cutoff( # 2019 Jan 08
day_max_temp,
day_min_temp,
threshold_temp,
cutoff_temp,
):
u"""Implement an intermediate cutoff.
"""
dd_threshold = gdd_single_sine_no_cutoff(day_max_temp, day_min_temp, threshold_temp)
dd_cutoff = gdd_single_sine_no_cutoff(day_max_temp, day_min_temp, cutoff_temp)
if None in [dd_threshold, dd_cutoff]:
result = None
else:
result = dd_threshold  2.0 * dd_cutoff
if result < ZERO:
result = ZERO
return result
def make_table(scale='F', species='Cydia pomonella'): # 2019 Jan 08
def f_to_c(f):
return int(round((f  32.0) * 5.0 / 9.0))
result = []
result.append(u'Growing DegreeDay (°%s) Table' % scale)
result.append(NULL)
range_max_temp = range(48, 120, 2)
range_min_temp = range(44, 92, 2)
if species in ['Cydia pomonella']:
common_name = 'Codling Moth'
threshold = 50
cutoff = 88
method = gdd_single_sine_horizontal_cutoff
elif species in ['Choristoneura rosaceana']:
common_name = 'Obliquebanded Leafroller'
threshold = 43
cutoff = 85
method = gdd_single_sine_vertical_cutoff
else:
raise NotImplementedError
if scale in ['F']:
pass
elif scale in ['C']:
range_max_temp = f_to_c(range_max_temp)
range_min_temp = f_to_c(range_min_temp)
threshold = f_to_c(threshold)
cutoff = f_to_c(cutoff)
else:
raise NotImplementedError
result.append(u'For %s (%s)' % (common_name, species))
result.append(u'Threshold = %3i°%s' % (threshold, scale))
result.append(u'Cutoff = %3i°%s' % (cutoff, scale))
result.append(u'Method: %s' % method.__name__)
result.append(NULL)
result.append(' \ Min Daily Temp')
cols = []
cols.append('Max')
for min_temp in range_min_temp:
cols.append(FMT % round(min_temp))
result.append(SPACE.join(cols))
for max_temp in range_max_temp:
cols = []
cols.append(FMT % round(max_temp))
for min_temp in range_min_temp:
dd_species = method(max_temp, min_temp, threshold, cutoff)
if dd_species is None:
cols.append('***')
else:
cols.append(FMT % round(dd_species))
result.append(SPACE.join(cols))
result.append(NULL)
return result
def main_line():
lines = make_table()
STDOUT.write('\n'.join(lines))
return
if __name__ == "__main__":
import doctest
doctest.testmod()
main_line()
# Fin
