6.6
"dd_table" Script
#!/usr/bin/python
# -*- coding: utf-8 -*-
# dd_table.py
# 2018 Feb 23 . ccr
u"""Print Growing Degree-Day Table.
This script generates a table of growing degree-days used for modeling
the emergence of larvae of Cydia pomonella (codling moth).
"""
# 2019 Mar 24 . ccr . Diagnose bad calculation.
# 2019 Jan 08 . ccr . Change function names. Provide
# . . gdd_single_sine_vertical_cutoff entry point
# . . among others.
# 2018 May 28 . ccr . Compare degree-days of cooling.
from __future__ import division
import sys
import codecs
import math
ZERO = 0
SPACE = ' '
NULL = ''
NUL = '\x00'
NA = -1
STDOUT = codecs.getwriter('utf-8')(sys.stdout)
FMT = '%3i'
def gdd_single_sine_with_theta( # 2019 Jan 08
day_max_temp,
day_min_temp,
threshold_temp,
):
u"""Return growing degree-days (GDD) above threshold, given daily max
and min temps.
This function is transcribed from a *perl* subroutine:
o Coop, Len. "Degree-Day 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 degree-days 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): 514-517.
I wouldn't say 1969 was early days for general-purpose 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 degree-days 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 real-world 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 degree-day 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 degree-days will (and should) outrun
the conventional degree-days.
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 24-hour (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 12-hour
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
try: # 2019 Mar 24
theta = math.atan2(d2, math.sqrt(day_diff_temp * day_diff_temp - d2 * d2))
except ValueError:
theta = math.pi
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 degree-days of cooling (DD).
This is a heating, ventilation, and air-conditioning (HVAC) rule of
thumb and is often held out as an adequate and easily grasped
substitute for calculating growing degree-days.
"""
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 Degree-Day (°%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
|