# Numerical Integration of Irregular Temperature Readings

I run a website that generates degree days, a specialist form of weather data used for calculations relating to building energy consumption. Without going into unnecessary detail, temperature is a function of time, and degree days are essentially the integral of that function.

At the moment our site calculates the data using an approximation method based on daily average, maximum, and minimum temperatures. But I'm working on improving that method by using finer-grained temperature measurements.

These finer-grained temperature measurements are taken throughout each day, and the recording interval can be anything from 1 minute to several hours - it depends on the weather station making the recordings. For most weather stations, the readings are fairly regular, but they're not completely regular. A weather station might typically record the temperature every half hour, but there will often be extra readings or missing readings, or readings taken at less regular intervals for certain periods.

Initially I've been using the trapezoidal method to numerically integrate the function of temperature against time. It's working pretty well, but I'm wondering if I might be able to improve it.

I'm not a mathematician, and my understanding of numerical integration is only very basic. I understand that Simpson's 1/3 rule and Simpson's 3/8 rule typically work better than the trapezoidal rule when numerically integrating mathematical functions. But real-world temperature readings don't follow an exact mathematical function. Also I understand that Simpson's rules require equal intervals, which my temperature readings don't consistently have.

I wonder if it might make sense to use Simpson's rules to integrate stretches of temperature readings that have 2 or more consecutive time intervals of equal length, and use the trapezoidal rule for stretches of irregular readings. But then I see here (a paper that I don't pretend to understand properly) that the trapezoidal rule can often work better than Simpson's rule for various classes of "rougher" functions. I would guess that outside air-temperature variation would be classed as pretty rough - the temperature jumps up and down throughout the day for all sorts of reasons.

I could probably come up with some way to estimate the effectiveness of various methods, but it's tricky because there's no "right answer" to compare figures against. So I'm trying to figure out what method would make most sense from a theoretical standpoint.

Do you think the trapezoidal rule is likely to be the best approach for me? Or are there other approaches that might make more sense?

• Can you give us order of magnitudes of the time between successive measurement, specially in the case where they are not equispaced ? The temperature have an obvious 24h cycle (and another 1 year cycle) which can pose problems in your integration if your measurement points are more than 12h appart... – Frédéric Grosshans Dec 3 '10 at 14:55
• OK reread the question, and you gave some answer in the text. My question will be : do you have to take account of duration longer than 6h without data ? – Frédéric Grosshans Dec 3 '10 at 14:57
• We're dealing with temperature data from thousands of weather stations, and they vary a lot, so it's difficult to generalize. If there are no readings for 24 hours we just estimate the degree days for that day by averaging the degree days from surrounding days. That's probably about the best we can do in that instance. Most stations have readings at least every 3 hours, but a minority of stations only measure temperatures between, say, 6 AM and 6 PM. Clearly not good for degree-day calculations because you miss the overnight temperature drop. We might just not use those stations at all. – MB01 Dec 3 '10 at 16:10

Since your temperature readings aren't equispaced, you can't directly apply Simpson's rule; the approach equivalent to this is to construct the parabola that interpolates three consecutive points (i.e., across two panels), and integrate that. The problem with this approach, of course, is that you need to have an odd number of data points (even number of panels) to do this.

You can use the trapezoidal rule, of course, across each panel, but a probably better idea might be to construct a cubic Hermite interpolant for each panel, and then integrate that. An obvious problem is that four conditions are needed to uniquely determine a cubic for each panel (two points and two derivative values), but estimates of derivative values can be constructed from the data such that the piecewise interpolant is locally monotonic; briefly, an piecewise interpolant is locally monotonic if there are no spurious inflection points or extrema within a panel. One approach to estimating derivative values for monotonic interpolation, due to Fritsch and Carlson, is implemented in FORTRAN as the pchip package, and in MATLAB as the function pchip. A more modern approach, and one which may be better is some situations, is due to Steffen. You may have to experiment which of the trapezoidal rule, Fritsch-Carlson, or Steffen would be best for integrating your data.

Let me detail the way one uses cubic Hermite interpolation for integrating data:

Each panel is bounded by two points, $(x_i,y_i)$ and $(x_{i+1},y_{i+1})$. The trapezoidal rule consists of constructing the line joining these two points (linear interpolation):

$$f_i(x)=\frac{x_{i+1}-x}{x_{i+1}-x_i}y_i+\frac{x-x_i}{x_{i+1}-x_i}y_{i+1}$$

and integrating that:

$$\int_{x_i}^{x_{i+1}}f_i(x)\mathrm dx=\frac{x_{i+1}-x_i}{2}(y_i+y_{i+1})$$

Integrating with cubic Hermite interpolation can be considered as a further "improvement" of the trapezoidal rule; briefly, in addition to points, one has derivative values $y_i^{\prime}$ and $y_{i+1}^{\prime}$, from which one constructs the cubic

$$g_i(x)=y_i+y_i^{\prime}(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3$$

where

$\displaystyle c_i=\frac{3\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-2y_i^{\prime}-y_{i+1}^{\prime}}{x_{i+1}-x_i}$ and $\displaystyle d_i=\frac{y_i^{\prime}+y_{i+1}^{\prime}-2\frac{y_{i+1}-y_i}{x_{i+1}-x_i}}{(x_{i+1}-x_i)^2}$.

Integrating $g_i(x)$ might look slightly complicated, however, there is a nice expression for the integral:

$$\int_{x_i}^{x_{i+1}}g_i(x)\mathrm dx=\frac{x_{i+1}-x_i}{6}\left(y_i+4g_i\left(\frac{x_{i+1}+x_i}{2}\right)+y_{i+1}\right)$$

whose verification I'll leave up to you.

I'll just add the note that in the case of equispaced data, integration with a piecewise cubic Hermite interpolant is equivalent to integration with the trapezoidal rule plus corrections at the beginning and end.

• Thank you J.M. - you've given me several things to look into and play about with. Much appreciated :) – MB01 Dec 3 '10 at 16:36
• The only warning would be the mantra "High order does not mean high accuracy". If you interpolate using higher order rules, that will improve accuracy if the real function looks like a polynomial and make it worse if it does not. This is not to detract from the excellent answer of J.M., just to say that in your case the trapezoidal rule may be at least good enough. – Ross Millikan Dec 4 '10 at 2:42
• @Ross is correct here @MB; like I said, you need to experiment and compare. Ensure that whatever result you get is physically justifiable! – J. M. is a poor mathematician Dec 4 '10 at 2:45
• Thank you both - that's very useful. I'm pretty sure the trapezoidal approach is going to be absolutely fine for stations that read the temperature every 5 minutes or so. It's the stations with many less readings (60 or 120+ minute intervals on average) that I wonder about. I've gathered some sample 5-minute data and my plan is to assume the trapezoidal result from that data as the "right answer". Then I'll strip out most of the readings, reducing it to 1-hourly or 2-hourly data, and try these alternative approaches on that lower-resolution data to see how the results compare. – MB01 Dec 4 '10 at 11:12