Guest post by Nick Stokes,
Every now and then, in climate blogging, one hears a refrain that the traditional min/max daily temperature can’t be used because it “violates Nyquist”. In particular, an engineer, William Ward, writes occasionally of this at WUWT; the latest is here, with an earlier version here. But there is more to it.
Naturally, the more samples you can get, the better. But there is a finite cost to sampling limitation; not a sudden failure because of “violation”. And when the data is being used to compile monthly averages, the notion promoted by William Ward that many samples per hour are needed, that cost is actually very small. Willis Eschenbach, in comments to that Nyquist post, showed that for several USCRN stations, there was little difference to even a daily average whether samples were every hour or every five minutes.
The underlying criticism is of the prevailing method of assessing temperature at locations by a combined average of Tmax and Tmin = (Tmax+Tmin)/2. I’ll call that the min/max method. That of course involves just two samples a day, but it actually isn’t a frequency sampling of the kind envisaged by Nyquist. The sampling isn’t periodic; in fact we don’t know exactly what times the readings correspond to. But more importantly, the samples are determined by value, which gives them a different kind of validity. Climate scientists didn’t invent the idea of summarising the day by the temperature range; it has been done for centuries, aided by the min/max thermometer. It has been the staple of newspaper and television reporting.
So in a way, fussing about regular sample rates of a few per day is theoretical only. The way it was done for centuries of records is not periodic sampling, and for modern technology, much greater sample rates are easily achieved. But there is some interesting theory.
In this post, I’d like to first talk about the notion of aliasing that underlies the Nyquist theory, and show how it could affect a monthly average. This is mainly an interaction of sub-daily periodicity with the diurnal cycle. Then I’ll follow Willis in seeing what the practical effect of limited sampling is for the Redding CA USCRN station. There isn’t much until you get down to just a few samples per day. But then I’d like to follow an idea for improvement, based on a study of that diurnal cycle. It involves the general idea of using anomalies (from the diurnal cycle) and is a good and verifiable demonstration of their utility. It also demonstrates that the “violation of Nyquist” is not irreparable.
Here is a linked table of contents:
- Aliasing and Nyquist
- USCRN Redding and monthly averaging
- Using anomalies to gain accuracy
- Conclusion
Aliasing and Nyquist
Various stroboscopic effects are familiar – this wiki article gives examples. The math comes from this. If you have a sinusoid frequency f Hz (sin(2π)) samples at s Hz, the samples are sin(2πfn/s), n=0,1,2… But this is indistinguishable from sin(2π(fn/s+m*n)) for any integerm (positive or negative), because you can add a multiple of 2π to the argument of sin without changing its value.
But sin(2π(fn/s+m*n)) = sin(2π(f+m*s)n/s) that is, the samples representing the sine also representing a sine to which any multiple of the sampling frequency s has been added, and you can’t distinguish between them. These are the aliases. But if s is small, the aliases all have higher frequency, so you can pick out the lowest frequency as the one you want.
This, though, fails if f>s/2, because then subtracting s from f gives a lower frequency, so you can’t use frequency to pick out the one you want. This is where the term aliasing is more commonly used, and s=2*f is referred to as the Nyquist limit.
I’d like to illuminate this math with a more intuitive example. Suppose you observe a running track, circle circumference 400 m, from a height, through a series of snapshots (samples) 10 sec apart. There is a runner who appears as a dot. He appears to advance 80 m in each frame. So you might assume that he is running at a steady 8 m/s.
But he could also be covering 480m, running a lap+80 between shots. Or 880m, or even covering 320 m the other way. Of course, you’d favour the initial interpretation, as the alternatives would be faster than anyone can run.
But what if you sample every 20 s. Then you’d see him cover 160 m. Or 240 m the other way, which is not quite so implausible. Or sample every 30 s. Then he would seem to progress 240m, but if running the other way, would only cover 160m. If you favour the slower speed, that is the interpretation you’d make. That is the aliasing problem.
The critical case is sampling every 25s. Then every frame seems to take him 200m, or halfway around. It’s 8 m/s, but could be either way. That is the Nyquist frequency (0.04 Hz), relative to the frequency 0.02Hz which goes with as speed of 8 m/s. Sampling at double the frequency.
But there is one other critical frequency – that 0.2 Hz, or sampling every 50s. Then the runner would appear not to move. The same is true for multiples of 50s.
Here is a diagram in which I show some paths consistent with the sampled data, over just one sample interval. The basic 8 m/s is shown in black, the next highest forward speed in green, and the slowest path the other way in red. Starting point is at the triangles, ending at the dots. I have spread the paths for clarity; there is really only one start and end point.
All this speculation about aliasing only matters when you want to make some quantitative statement that depends on what he was doing between samples. You might, for example, want to calculate his long term average location. Now all those sampling regimes will give you the correct answer, track centre, except the last where sampling was at lap frequency.
Now coming back to our temperature problem, the reference to exact periodic processes (sinusoids or lapping) relates to a Fourier decomposition of the temperature series. And the quantitative step is the inferring of a monthly average, which can be regarded as long term relative to the dominant Fourier modes, which are harmonics of diurnal. So that is how aliasing contributes error. It comes when one of those harmonics matches the sample rate.
USCRN Redding and monthly averaging
Willis linked to this NOAA site (still working) as a source of USCRN 5 minute AWS temperature data. Following him, I downloaded data for Redding, California. I took just the years 2010 to present, since the files are large (13Mb per station per year) and I thought the earlier years might have more missing data. Those years were mostly gap-free, except for the last half of 2018, which I generally discarded.
Here is a table for the months of May. The rows are for sampling frequencies of 288, 24, 12, 4, 2, and 1 per day. The first row shows the actual mean temperature averaged 288 times per day over the month. The other rows show the discrepancy for the lower rate of sampling, for each year.
| Per hour | 2010 | 2011 | 2012 | 2013 | 2014 | 2015 | 2016 | 2017 | 2018 |
| 1/12 | 13.611 | 14.143 | 18.099 | 18.59 | 19.195 | 18.076 | 17.734 | 19.18 | 18.676 |
| 1 | -0.012 | 0.007 | -0.02 | -0.002 | -0.021 | -0.014 | -0.007 | 0.002 | 0.005 |
| 2 | -0.004 | 0.013 | -0.05 | -0.024 | -0.032 | -0.013 | -0.037 | 0.011 | -0.035 |
| 6 | -0.111 | -0.03 | -0.195 | -0.225 | -0.161 | -0.279 | -0.141 | -0.183 | -0.146 |
| 12 | 0.762 | 0.794 | 0.749 | 0.772 | 0.842 | 0.758 | 0.811 | 1.022 | 0.983 |
| 24 | -2.637 | -2.704 | -4.39 | -3.652 | -4.588 | -4.376 | -3.982 | -4.296 | -3.718 |
As Willis noted, the discrepancy for sampling every hour is small, suggesting that very high sample rates aren’t needed, even though they are said to “violate Nyquist”. But they get up towards a degree for sampling twice a day, and once a day is quite bad. I’ll show a plot:
The interesting thing to note is that the discrepancies are reasonably constant, year to year. This is true for all months. In the next section I’ll show how to calculate that constant, which comes from the common diurnal pattern.
Using anomalies to gain accuracy
I talk a lot about anomalies in averaging temperature globally. But there is a general principle that it uses. If you have a variable T that you are trying to average, or integrate, you can split it:
T = E + A
where E is some kind of expected value, and A is the difference (or residual, or anomaly). Now if you do the same linear operation on E and A, there is nothing gained. But it may be possible to do something more accurate on E. And A should be smaller, already reducing the error, but more importantly, it should be more homogeneous. So if the operation involves sampling, as averaging does, then getting the sample right is far less critical.
With global temperature average, E is the set of averages over a base period, and the treatment is to simply omit it, and use the anomaly average instead. For this monthly average task, however, E can actually be averaged. The right choice is some estimate of the diurnal cycle. What helps is that it is just one day of numbers (for each month), rather than a month. So it isn’t too bad to get 288 values for that day – ie use high resolution, while using lower resolution for the anomalies A, which are new data for each day.
But it isn’t that important to get E extremely accurate. The idea of subtracting E from T is to remove the daily cycle component that reacts most strongly with the sampling frequency. If you remove only most of it, that is still a big gain. My preference here is to use the first few harmonics of the Fourier series approximation of the daily cycle, worked out at hourly frequency. The range 0-4 day-1 can do it.
The point is that we know exactly what the averages of the harmonics should be. They are zero, except for the constant. And we also know what the sampled value should be. Again, it is zero, except where the frequency is a multiple of the sampling frequency, when it is just the initial value. This is just the Fourier series coefficient of the cos term.
Here is are the corresponding discrepancies of the May averages for different sampling rates, to compare with the table above. The numbers for 2 hour sampling have not changed. The reason is that the error there would have been in the 8th harmonic, and I only resolved the diurnal frequency up to 4.
| Per hour | 2010 | 2011 | 2012 | 2013 | 2014 | 2015 | 2016 | 2017 | 2018 |
| 1/12 | -0.012 | 0.007 | -0.02 | -0.002 | -0.021 | -0.014 | -0.007 | 0.002 | 0.005 |
| 2 | -0.004 | 0.013 | -0.05 | -0.024 | -0.032 | -0.013 | -0.037 | 0.011 | -0.035 |
| 6 | 0.014 | 0.095 | -0.07 | -0.1 | -0.036 | -0.154 | -0.016 | -0.058 | -0.021 |
| 12 | -0.062 | -0.029 | -0.075 | -0.051 | 0.019 | -0.066 | -0.012 | 0.199 | 0.16 |
| 24 | 1.088 | 1.021 | -0.665 | 0.073 | -0.864 | -0.651 | -0.258 | -0.571 | 0.007 |
And here is the comparison graph. It shows the uncorrected discrepancies with triangles, and the diurnally corrected with circles. I haven’t shown the one sample/day, because the scale required makes the other numbers hard to see. But you can see from the table that with only one sample/day, it is still accurate within a degree or so with diurnal correction. I have only shown May results, but other months are similar.
Conclusion
Sparse sampling (eg 2/day) does create aliasing to zero frequency, which does affect accuracy of monthly averaging. You could attribute this to Nyquist, although some would see it as just a poorly resolved integral. But the situation can be repaired without resort to high frequency sampling. The reason is that most of the error arises from trying to sample the repeated diurnal pattern. In this analysis I estimated that just from Fourier series of hourly readings from a set of base years. If you subtract a few harmonics of the diurnal, you get much improved accuracy for sparse sampling of each extra year, at the cost of just hourly sampling of a reference set.
Note that this is true for sampling at prescribed times. Min/max sampling is something else.
Finally somebody objects to the nonsense of the author who decided to bark at the Nyquist tree.
“Min/max sampling is something else. – Nick Stokes”
Nick,
Thank for the detailed presentation, it makes good sense to me.
I wish you had (Or in a future post), looked into Min & Max as a special case of “sampling” though.
It is interesting to note that with only a Min & Max sample, your “runner on the track” could only ever be found where she actually was, at the start or the halfway mark* of the circuit. This aperiodic cycle “selection” could at worst, underestimate the “squareness” of the diurnal wave shape if “he” happened to be resting – hare like – at those locations during the race! The time between points is always 12h though. Sot that lengthening the time spent at either extreme also reduces the time between them.
cheers,
Scott
*Top and bottom, peak and trough of the cycle
**See: Aesop’s Fable, The Tortoise and the Hare
Grrrrrr! Typos! I only ever see them after posting! ‘My kingdom’ for an editing function! 😉
I guess the problem I’m having is the simple fact that a day with T-min of 40f and T-max of 80f will average out to be exactly the same as a day with T-min of 55f and T-max of 65f. The actual experience at that location would be very different on those two days, but would not end up being represented in the final average. Is the actual daily temperature range in a given location considered insignificant?
Yeah I agree,
I wasn’t talking about Tmean = (max+min)/2 , just the comparison between 2 random samples and the selection of the 2 extremes max and min!
Nick, you wrote above –
“Most of the error in monthly average results from the interaction of the sampling frequency with harmonics of the diurnal frequency. So if you subtract even a rough version of the latter, you get rid of much of the error.”
This is the assertion that worries me, one that I find counter-intuitive. You should not be able to improve accuracy by mathematics alone. I think you are improving an illusion of accuracy, but that has no significant value.
Indeed, it raises another point of whether it is even valid to use Fourier Transform concepts and Nyquist concepts. There are classes of stochastic data where its is inappropriate. What if we are looking a data filling the main requirements of a bounded random walk? These daily temperatures are highly autocorrelated with next day.
Geoff.
Geoff,
“You should not be able to improve accuracy by mathematics alone.”
The situation with say 2/day sampling is that you have 60 readings to estimate a month average. That should be enough, but there is a problem that a repetitive diurnal variation can interact with the sampling rate to produce an error. There is no reason why maths can’t identify and repair that error. It is the same every time – not particular to the latest month’s data.
Nick,
There is so much wrong with your answer that I am going to sit on it to cogitate about how to express how wrong your assumptions are.
First, though, are these historic temperature numbers in a class of numbers that allow valid Fourier analysis? I doubt they are. Geoff.
Geoff,
Actually, I’ve realised I don’t even need Fourier analysis for the improvement. But anyway it was applied to the estimated periodic diurnal variation, so yes, it is valid. I’ll include the non-Fourier version in my response to WW.
While math alone cannot consistently improve the accuracy of individual measurements, it’s not powerless in correcting for well-understood sources of systematic BIAS in the AVERAGES of data. The asymmetry of the AVERAGE diurnal wave-form, an expression of phase-locked harmonics, is a well-documented, PERSISTENT source of discrepancy between the daily mid-range value and the true mean. Since that asymmetry is produced by factors driven by astronomically repeating patterns of insolation, present-day determinations of the resulting discrepancy provide valid corrections for the station-dependent bias in historical monthly averages.
Questioning the validity of universally applicable DSP methods without any analytic understanding of them is an empty exercise. Especially so when highly autocorrelated, quasi-periodic data in all scientific fields of study are duck soup for proper DSP techniques.
Nick Stokes said:
“Naturally, the more samples you can get, the better. But there is a finite cost to sampling limitation; not a sudden failure because of “violation”.”
Thanks Nick. I’d been trying to get my head around this for a while, and I think I’ve now got a better idea of how much this matters and why.
While Nyquist only was conceived in terms of a repetitive signal the same error can be seen to be produced spatially. The number of stations needed for global warming suddenly multiplies by a thousand or more once you accept the distribution of CO2 from fossil fuel and the temperatures is no longer even.
Also the effect of fossil fuel is considered to be a warming one so we are actually talking energy not temperature so Hoyt Clagwell has a point in that a reduction in the lower temperature looks like an increase in the apparent fossil fuel induced energy.
Hi Nick,
Thanks for the really interesting post. You have done some interesting work here and I’m glad we have the opportunity to explore this further. I have some comments and questions. But my goal is to get to your exercise in correcting for the aliasing – that is the important thing here.
First, can you explain the first column of the charts for your Redding data? Should the heading be “Hours between samples”? If so, then I think I understand what you are doing, but “per hour” has me confused.
Nick said: “Willis Eschenbach, in comments to that Nyquist post, showed that for several USCRN stations, there was little difference to even a daily average whether samples were every hour or every five minutes.” And: “As Willis noted, the discrepancy for sampling every hour is small, suggesting that very high sample rates aren’t needed, even though they are said to “violate Nyquist”. But they get up towards a degree for sampling twice a day, and once a day is quite bad.”
My reply: Willis and I are actually in reasonably good agreement at this point. To summarize: He has shown that with averaging USCRN data, 24-samples/day seems to be the point where the error stops reducing “significantly” as sample rate is increased. With the limited study done I have no quarrel with this. However, Engineers design systems not for the average condition, but for the extreme conditions. I showed the example using USCRN data for Cordova AK, where another +/-0.1C error could be reduced by going to 288-samples/day. 288 is what NOAA uses in USCRN so it is reasonable to select this number, but further research could justify increasing or decreasing this number. We don’t know what kind of spectral profiles we would see if we sampled at other locations. A few have mentioned the data available from Australia as being sampled every 1 second. Perhaps I’ll have a look at that data someday. Additionally, from an engineering perspective we want additional margin. A higher sample rate makes the anti-aliasing filter easy to implement. And the cost of sampling higher, both for the equipment and to manage the data is arguably insignificant. If sampled properly, it is easy to reduce data and filter out frequencies digitally, with no negative consequences.
Nick said: “That of course involves just two samples a day, but it actually isn’t a frequency sampling of the kind envisaged by Nyquist. The sampling isn’t periodic; in fact, we don’t know exactly what times the readings correspond to. But more importantly, the samples are determined by value, which gives them a different kind of validity. Climate scientists didn’t invent the idea of summarising the day by the temperature range; it has been done for centuries, aided by the min/max thermometer. It has been the staple of newspaper and television reporting.
My reply: This has been a sticking point in the conversations. Nyquist tells us how to make sure any digital domain representation of an analog domain signal maintains its relationship with the analog signal. I think we all agree that the goal of working with a set of discrete measurements of a continuous signal is to get a result that is for the analog signal that happened in the real world. It doesn’t matter how we come by the samples. Any discrete value measurement of an analog signal is a sample. If the samples do not meet the requirements to reconstruct the original, then Nyquist cannot be complied with – said another way it is violated. Max and min are periodic: These samples happen 2x/day every day. They just happen with error in the periodicity.
Nick said: “The interesting thing to note is that the discrepancies are reasonably constant, year to year. This is true for all months. In the next section I’ll show how to calculate that constant, which comes from the common diurnal pattern.”
My reply: Nick, I studied 26 stations for trends. I agree with you from this limited study, that the offset from the high sample rate reference tends to be constant from year-to-year. However, each station looks distinctly different from the others. Some stations show very little offset but large trend errors (Ex: Baker CA). Some show large offsets (1.5C!) but very small trend errors. (Ex: Fallbrook CA) Some stations have positive offsets and some negative. And some have offsets that change sign every few years (Ex: Montrose CO). Studying many stations may be beneficial.
https://imgur.com/IC7239t
https://imgur.com/cqCCzC1
https://imgur.com/SaGIgKL
https://imgur.com/xA4hGSZ
Regarding your section about “Using Anomalies to Gain Accuracy”: Can you explain in more detail what you actually did? I see that your second table shows better results, but what did you do exactly?
Nick said: “But it isn’t that important to get E extremely accurate. The idea of subtracting E from T is to remove the daily cycle component that reacts most strongly with the sampling frequency.
My reply: How do you determine the daily cycle component from an aliased signal?
Fundamentally speaking, aliasing is the transformation of high frequency information into false low frequencies that were not present in the original signal. In the example of 2-samples/day, frequency content at 2-cycles/day shows up as energy at 0-cycles/day – the “zero-frequency” or “DC”. We are dealing with samples – so it is not possible to distinguish the original DC energy from the aliased DC energy. The energy combines in your samples. Aliasing is generally considered irreversible and not repairable – unless you have some very specific information and know what to remove. I think you are saying that this is what you are doing – but I’d like to know the details. What is your source of data for this information, regarding Redding what are the numbers and what is the operation you perform? Are you using information from 288-samples/day? If so, I think this is ironic because if we sample correctly, we don’t need to correct. Can you correct a record where you don’t have high sample rate data to support you? I suggest that you take max and min values from the 288, discard the timing information and then do your correction method to come up with a signal that is close to that provided by 288-samples/day. We can use the 288 signal as reference to check your work, but you can’t have that data to get the fix. For this to be meaningful, it has to be usable on the historical record that doesn’t benefit from an alternate record. Also, I recommend you try this on Blackville SC, Monahans TX, or Spokane WA in USCRN (I can provide the files). They present the largest trend errors in my study. With Montrose CO, the offset sign changes twice over a 12 year period. Would your method hold up to this condition?
The trend errors I show are of similar magnitude to the claimed trends we are experiencing. If we had a properly sampled record, would those warming trends be reduced or increased? We really don’t know. So far, we didn’t get to any critical review of the 26 station trend errors presented.
In your conclusion you say: “Sparse sampling (eg 2/day) does create aliasing to zero frequency, which does affect accuracy of monthly averaging. You could attribute this to Nyquist, although some would see it as just a poorly resolved integral. But the situation can be repaired without resort to high frequency sampling.”
My response: It would be great if the “can be repaired without resort to high frequency sampling” is true. I’m not ready to agree to that yet Nick. I really don’t know what you did. For something this important – because it would be a big thing – I’d like to see the solution work on some tougher test cases with all of the work shown. If it works would you be willing to process some annual records (max/min) for one of the datasets to see how the yearly anomaly is affected?
Thanks, William
It’s been a long day here – Australia day stuff. I’ll respond tomorrow
REgards,
Nick
Hey Nick,
I was wondering where you were from. Are you living in Australia? I saw that occasionally you spelled a word “British style”, with an S, where Americans use a Z. Ex: organization vs organisation. I’m in the US (Atlanta). I hope you had a good Australia Day!
Take your time on the reply and let me know if it isn’t clear what I’m asking. If you can provide me a quick tutorial about how to apply your correction, I can do the work on my end to validate with test cases if you prefer.
Dr. Stokes is indeed Australian. He worked for over 35 years as research scientist in the Division of Mathematics and Statistics at the CSIRO.
Hello Philip,
Thanks for the reply and the information. It sounds like a very distinguished career. I found a video about an award that Nick and his team won for the development of “Fastflo” computational fluid dynamics software, in 1995. Very impressive stuff. I’m sure there is much more.
I also looked up Canberra, where the CSIRO is located. It appears to be an absolutely beautiful place.
William,
CSIRO is spread out; I’m in Melbourne (though I started in Canberra). The Australia day event was the open day at our historic Government House, where they put on quite a show. I try to round up as many grandchildren as I can.
Haven’t quite got the new, Fourier-free code working yet. It will be easier to replicate when I have.
Thanks for the reply Nick. I was just reading an article about how Melbourne is closing in on Sydney as Australia’s largest city. Also mentioned was a 356-meter tall building that would be Australia’s tallest. I guess Australia has a lot of beautiful places as the images I have seen of the Melbourne skyline with its nearby parks looks very inviting. I know Australia has a lot of unforgiving land as well.
I’m glad you have a large family to share the holiday with.
I’ll look forward to your reply once you get your code working.
William
“If you can provide me a quick tutorial about how to apply your correction, I can do the work on my end to validate with test cases if you prefer.”
Here is the revised method, which is more straightforward, doesn’t use FFT, and seems to give even better results. Firstly, the basic algebra again:
We have a high res linear operator L, which here is monthly averaging with hourly sampling, and a low-res, which is monthly averaging with 12,6,2 or 1 samples per day. I split T into E and A. E is an estimate of average diurnal for the month.
So the improvement on L1(T) is
L1(A) + L(E) = L1(E+A) + L(E) – L1(E)
The point of writing it this way is that L1(E+A)=L1(T) is just the low res calc, and L(E) – L1(E) just depends on the 12 sets of hourly diurnal estimates – a 24 x 12 matrix.
I start with an array T(hour, day, month, year). That is 24 x 31 x 12 x 9 (9 years 2010-2018).
Then
1) Get E(hour, month) by averaging over 31 days and 5 years (2010-2014). You could use fewer years as reference. The result E is 24 x 12
2) Get L1(T) by decimating the hour part of T – eg for 2/day by taking hours 1 and 13. That gives a T1 array 2 X 31 x 12 x 9. Then average over the decimated hours and days, to get a 12 x 9 set of monthly estimates
3)E is 24 x 12. Get L(E) by averaging over all hours, and L1(E) by averaging over the decimated set. Then form L(E) – L1(E) – just 1 number for each month.
4) So for the corrected set, just add L(E) – L1(E) to each month of L1(T).
Here is my current R code
#source(“c:/mine/starter.r”); rung(“C:/mine/blog/data/uscrn”,”approx.r”)
# reads 5 min USCRN dada from ftp://ftp.ncdc.noaa.gov/pub/data/uscrn/products/subhourly01
s=”ftp://ftp.ncdc.noaa.gov/pub/data/uscrn/products/subhourly01/2018/CRNS0101-05-2018-CA_Redding_12_WNW.txt”
ave=function(x,k){ # averaging multidim array; averages dims where k>0
d=dim(x);
for(i in 1:length(d)){
j=k[i]
if(j>0){
x=apply(x,-i,mean,na.rm=T);
d[i]=1; dim(x)=d
}
}
x
}
if(0){ # read month files from web and extract – run once only
data5min=NULL
for(i in 2010:2019){
f=gsub(“2018”,i,s)
b=readLines(f)
nb=length(b)
v=matrix(“”,nb,3)
v[,1]=substr(b,21,29)
v[,2]=substr(b,30,34)
v[,3]=substr(b,58,65)
u=as.numeric(v)
dim(u)=dim(v)
data5min=rbind(data5min,u)
}
save(data5min,file=”Redding.rda”)
}###
if(1){#sort 5min data into array by month
x=data5min[,1]
j=which(diff(x)>50)
N=288*31
array5min=daysinmth=NULL
for(i in 2:length(j)-1){
ij=(j[i]+1):j[i+1]
daysinmth=c(daysinmth,length(ij))
array5min=cbind(array5min,pad(data5min[ij,3],N))
}
array5min[array5min< -99]=NA
dim(array5min)=c(288,31,12,9)# 5min,day,mth,yr
hrs=seq(1,288,12)
arrayhrs=array5min[hrs,,,] # reduce from 5min to hourly
hires=ave(arrayhrs,c(1,1,0,0))
diurnalhr=ave(arrayhrs[,,,1:5],c(0,1,0,1)) # hr x mth
}###
L=function(x){ # Hi res averaging (hourly)
ave(x,c(1,1,0,0))
}
L1=function(x,n=2){ # Low res averaging
d=dim(x)
i=seq(1,24,24/n)
d[1]=length(i)
str(x)
x=array(x[i,,,],d)
ave(x,c(1,1,0,0))
}
if(1){ #Calculates an array w of results
w=array(0,c(12,9,6))
iw=c(24,6,2,1)
w[,,1]=hires # The top row of hires
for(i in 1:4){ # for each low res, the difference from hires
approx=w[,,i+1]=(L1(arrayhrs,iw[i]) + c(L(diurnalhr)-L1(diurnalhr,iw[i])))-hires
}
w=round(w,3)
}
Now that is putting your money where your mouth is.
Thank you for the algorithm and code Nick.
WW,
You are into topics so seldom ventilated that terminology is ragged and even individual to the researcher.
Example, relevant to Australian MMTS equipment, synonym AWS, synonym PRT, approximately, mostly post 1 November 1996, the official date when this equipment was to become formally primary, but start dates before and after happen. Thus we have –
1 second samples. There is a lack of public documentation, but I take this to be the duration of the instrument’s data logger window, opened from time to time.
1 minute samples. Can be either a number of 1 second samples added, averaged, then reported each minute.
10 minute sampling. We have this is Aust. It seems to be recorded in 3 parts, the highest 1 second sample logged over 10 minutes, the lowest and the most recent. (The latter could be for some form of time marking, I do not know).
Daily sampling. Current practice is for the Tmax to be the highest 1 second sample in the midnight-midnight duration of a calendar day. The Tmin is the lowest such sample. However, there is uncertainty about this; it does not apply to the Liquid-in-glass period of thermometry. Importantly, the time of day when Tmax and Tmin were reached does not seem to be publicized routinely, though it might be archived. I do not know precisely what is archived.
So, you can see that we in Australia are essentially based on 1 second samples, using BOM designed PRT that has thermal inertia designed to match LIG response times. It follows that there is confusion about what “sampling” means. Mentally, I use the minimum time the data logger is acquiring voltage from the PRT. Conversationally, I am lax as the words above have different meaning for “sampling” This, of course, impacts in some ways on the (again loose) term “sampling frequency” which commonly seems to be a period longer than a second, over which temperatures are observed, as in “daily sampling” which terminology might vary from country to country.
What we export for the global data manufacturers like Hadley, NOAA, Best, etc is anybody’s guess. There are international standards available, but I do not know how well they are followed, nor how the job was done in the early, formational years. Also, we in Australia seem to have no control over what downstream adjusters do with our data. There is some clarification overall at
http://www.bom.gov.au/climate/cdo/about/about-directory.shtml
http://www.bom.gov.au/climate/cdo/about/about-airtemp-data.shtml
Currently, a broad brush look at Australian T data break into these broad parts –
pre-1910. Data judged unreliable for routine use by BOM.
1910-1992 approx. Automatic Weather Station commenced in remote places about 1992, followed by a rush in the late 1990s. The actual station dates are in
http://www.bom.gov.au/climate/data/stations/about-weather-station-data.shtml
1910- 1Sept 1972. Metrication, from degrees F to degrees C. This is the formal date. Some stations might have different dates in practice. Within this period there were some directions to read thermometers to only the nearest degree (no fractions).
1860s to now. A variety of screens was used around thermometers. They require corrections to be compatible. The record is scattered and a single national correction has not been available.
In summary, the plot thickens once we start on a concept of analysing sampling frequency at Australian stations. You simply cannot take 100 year slabs of data as being capable of mathematical analysis without either prior selection of sets of compatible stations, or adjustment of some sets to be compatible with others. Once adjustment commences, all Hell breaks loose.
Sorry to be so negative, but the errors inherent in the data are most often larger than the effect being sought (like how much global warming there was) and this poor signal:noise ratio dominates most forms of mathematical analysis. There is no blame game here. The old systems were designed for purposes different to those of today. They did their jobs, kudos to those involved, but they are simply unfit for purpose when that purpose is smaller in magnitude than the error. It is NOT OK to claim they are all we have, then proceed. It is past time to officially declare “UNFIT FOR PURPOSE” before say year 2000 and close the books on many future exercises in futility. Geoff.
Hello Geoff,
That is some really good information to consider about Australia. I looked at both links you provide from the BOM. I tried to find and download 1 minute data, but the system was temporarily unavailable. I’ll try again tomorrow. Meanwhile, if you can recommend some cities/stations with good 1-minute data I’d like to take a look at it and try to do some similar analysis to what I did for USCRN.
I agree that more information would be very helpful to understand the system – this is similar to the US. Others have provided me the links to the ASOS user guides but there are still many questions an engineer would ask that are not answered in the guides.
I have been working with the assumption that the datasets (GISS, HADCRUT) use the max/min method. Input from Nick and others seems to confirm that. So even though USCRN uses 20-sec samples (averaged to 5-min) and Australia has 1-sec samples, it doesn’t appear that this data is what goes into the datasets. So we have better data available, but it is not used for the global average calculations or trends. As you alluded to, it isn’t good to mix the new higher quality data with the older data, and it seems they don’t do this. It would be good to have parallel analysis using both methods for the years that the new method data is available.
I can give you my opinion about what 1-second samples means based upon my experience with Analog-to-Digital Converters (ADCs). First, these converter chips are designed for their intended applications. ADCs designed to sample communications signals at tens of megaHertz to tens or hundreds of gigaHertz will sample at multiples of those speeds. Instrumentation and process control need to sample signals in the kiloHertz or low megaHertz range (or slower). The “window” as you say, is open for a very small amount of time after the clock pulse arrives. You don’t want the signal to change much before you measure it. I think a good analogy is camera shutter speed. Faster shutter speeds are needed if the object you are trying to photograph is moving. If the shutter is fast enough compared to the motion of the object then the image on film is sharp. If the shutter is not fast enough then the image appears to blur and shows evidence of motion. So 1-second samples means that we take a picture every 1 second but the shutter speed is very fast compared to 1-second. The 1-sec data can be averaged to a slower rate if desired. Max and min can be determined from this 1-sec data. So a 1-sec sample rate means there are 60 samples in 1 minute and 3600 samples in 1 hour and 86,400 samples in each day. If a good quality ADC is used then we don’t need to concern ourselves with the window opening time. The converter takes care to get an accurate measurement before the signal can change value significantly.
I agree with you, for the reason of sampling and a dozen other reasons that the max/min record is not adequate to support the claims of imminent climate doom.
Hi WW,
I did some years of aeronautical engineering study at RAAF Academy before a car crash shortened my studies to a basic science degree because I could not longer fly. In working life I did a lot of both engineering and science as well as political motivation theory and the madness of crowds. So I guess you and I were destined to have a wavelength in harmony.
Bought my first ADC IN 1970, an expensive 8K frequency job in a nuclear radiation counter.
Thank you for your reminders that engineering success strongly depends on delivering the goods. Too much climate science these days is minor variations on set themes and too much of it fails to properly address error analysis. Best regards. Geoff.
Hello Geoff,
I’m sorry to hear about the car crash. Sometimes life gets redirected and who knows whether the original or new path is best. We can’t run 2 lives in parallel and choose the best one at the end. Hopefully, the slight redirection put you on the optimal path.
The radiation counter sounds like a an expensive toy – or perhaps this was a tool for work. We have a lot of amazing technology today, but also a lot of really good things since the 60’s and 70’s. Today we have many incremental improvements that rest on the innovation of that earlier work.
You mentioned “madness of crowds”. I can’t help but think of that phrase every time I watch the news about current events. An alien mind virus seems to have infected half of the population. While I didn’t read it, I’m aware of this publication (by Charles Mackay) and I assume you are familiar with it:
Extraordinary Popular Delusions and the Madness of Crowds:
https://en.wikipedia.org/wiki/Extraordinary_Popular_Delusions_and_the_Madness_of_Crowds
Some good quotes from it:
“Men, it has been well said, think in herds; it will be seen that they go mad in herds, while they only recover their senses slowly, and one by one.”
“We go out of our course to make ourselves uncomfortable; the cup of life is not bitter enough to our palate, and we distill superfluous poison to put into it, or conjure up hideous things to frighten ourselves at, which would never exist if we did not make them.”
“We find that whole communities suddenly fix their minds upon one object, and go mad in its pursuit; that millions of people become simultaneously impressed with one delusion, and run after it, till their attention is caught by some new folly more captivating than the first.”
Since you have studied the subject, what do you think of Mackay’s study? Is there a better reference recommended as a primer?
Thanks for the exchange Geoff and for your affirmations about the benefits engineering discipline could bring to climate science. I will look forward to future discussions with you for sure.
All the Best,
WW
The quaint notion that JUST ANY daily discrete value of an analog signal constitutes the sampling specified by the Sampling Theorem is readily dispelled even by Wiki-expertise:
https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem)
which patently requires an UNEQUIVOCALLY FIXED sampling rate fs of a Dirac comb. The longer this basic misconception is vehemently defended, the lesser the evidence of analytic competence.
Thank you 1sky1,
You argue against my comments while making my points for me. Quite a “quaint” talent you have.
Yes, perfect reconstruction requires perfect quantization and perfect (“UNEQUIVOCALLY FIXED”) periodic sampling. Any deviation from these perfections translates to imperfections or error in the reconstruct-ability and therefore imperfections or error in the sampling.
Nyquist is the “bridge” or transform between analog and digital domains. If you are working with samples that don’t comply with Nyquist then you violate Nyquist. Your samples don’t accurately represent the analog signal. So math operations on those samples poorly apply to the original analog signal. Isn’t that the entire point of working with samples?
Competence is not measured by reading Wikipedia, but decades of real world hands-on experience and billions of channels of converters deployed in the real world.
The astonishing blunder of mistaking metrics of the analog signal (such as Tmax and Tmin) as a product of “imperfections or error in the [discrete] sampling,” is now given the rhetorical cover of “reconstruct-ability.” That there’s absolutely nothing to reconstruct when those metrics are obtained in practice SOLELY from CONTINUOUSLY sensing thermometers math totally escapes the marketing mind that brags about “hands-on experience and billions of channels of converters deployed.”
The claim that “math operations on those samples poorly apply to the original analog signal” gives ample proof of chronic inability to recognize the fundamentals of an intrinsically different problem. The registering of extrema is a signal-dependent operation, independent of any predetermined “sampling” scheme. In that context, “violating Nyquist” is a total red herring.
1sky1 is quite correct. The Nyquist argument as applied to two different numbers (not samples – since we don’t have associated times) is wrongly offered as a cause of error in the calculation (Tmax+Tmin)/2 as a substitute for mean. The average is just an error that is clearly WRONG and does not need to violate any fancy law to disparage it.
1sky1 suggested January 23, 2019 at 4:02 pm on William’s thread an empirically- determined equation for the mean: (1 – eta)Tmin + eta Tmax, which is a constructive suggestion. Eta was not determined.
On William’s thread January 26, 2019 at 7:53 pm, not giving up on sampling (using a non-uniform approach): http://electronotes.netfirms.com/AN356.pdf , I suggested how we might put in some typical values of times for Tmax and Tmin (best approximations from a dense uniform time grid) and arrive at a firsts (theoretical) approximation for sky’s eta.
I was not optimistic about getting the work done, but was surprised to find that I HAD already (back in 2003) written code (Matlab/Octave) that I was dreading even thinking about.
I am almost certainly NOT going to finish this project (right or wrong) soon enough that
this thread will still be active. I will write up the results as one of my Electronotes Webnotes for anyone caring to check. Or email me: hutchins@ece.cornell.edu.
– Bernie
Bernie,
I looked at your note. It seems to me that if you restrict the problem to sampling a finite number of times per day, same times every day, then it is possible:
1. You can regard it as the sum of 1/day sampling at the various times. eg if you sample at 7am and 4pm, then it is the sum of 7am sampling and 4pm sampling.
2. So only harmonics of the diurnal can give DC components.
3. So if you have the Fourier coefficients for the diurnal, then the DC component is just the sum of the coefficients times the corresponding sampled trig functions
eg if you sample at t1 and t2 times of day, then the component of the Fourier series b2*cos(w*t), w=diurnal angular frequency, would yield DC b2*cos(w*t1)+b2*cos(w*t2)
etc
But I don’t think the eta approach is the right idea. The dependence of min/max on time of observation dominates, and you can’t get an answer without bringing that in. In sampling terms, if TOBS is say 5pm,for a reasonable number of days the time of max will swing from say 3pm of that day to 5pm the previous day, if that was warmer. That completely messes up sampling.
Nick – Thanks for looking at my app note : http://electronotes.netfirms.com/AN356.pdf
Let me go over this again from the perspective of someone such as yourself who has actually read my app note. [More to the point, I am trying to explain to myself as I go along!] The note relates a “sampling cell”, which I call SAMCELL, to a matrix M which describes exactly how an original spectrum is messed up by irregular sampling. The original spectrum is described by a corresponding “spectral cell” which I call SPECCELL. Both “cells” are the same length, L, as periodic sequences. SAMCELL is NOT a time signal, but rather a vector of placeholders for potential samples (1 if sample is taken, 0 otherwise), and SPECCELL is NOT a frequency spectrum but rather denotes equal-width segments on the interval of 0 to half the sampling frequency which may or may not be occupied (1 if occupied, 0 otherwise).
We begin with the notion that the daily 24-hour temperature cycle would approximate (perhaps not too well) a sinusoidal-like cool-down-warm-up cycle. We consider taking 24 hourly measurements. Nick has reasonably suggested that a min might often occur at the 7 AM hour while a max might often occur at 4 PM, and we take samples to be at those times, wishing ourselves good luck. Thus SAMCELL becomes, starting midnight: [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0] – very sparse.
Now – what is SPECCELL? SPECCELL can have at most 2 non-zero values, else we violate the general form of the sampling theorem (bunched sampling). SPECCELL indexes frequencies 0, 1/24, 2/24, 3/24 . . . 23/24 (as with a DFT). Thus we might need to filter or otherwise be assured of a greatly reduced bandwidth. For the present problem, however, we are not looking to decimate any known full time-sequence, but rather see HOW the spectrum of a highly bandlimited signal could be pieced back together. Thus we ASSUME that we are interested in reconstructing a spectrum already known or forced to be bandlimited as SPECCELL = [1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]. That is we are looking for a spectrum consisting of only a DC term and a sinusoidal of frequency 1/24 (one per day). This IS an artificial signal for which the mean IS the middle value between extremes, for which these extremes must be 12 hours apart. But we chose them to be at 7 AM and 4 PM – so to continue.
Here however we want to start with input values we have or can arguably assumes: Tmax and Tmin (given traditional measurements); the 7 AM and 4 PM (assumed times – SAMCELL above); and the supportable bandwidth (SPECCELL above). We would form the stand-in signal as Tmin at 4 AM and Tmax at 7 PM) and compute its spectrum, likely trimming the FFT to the first two values.
Now the heart of the development is to use the M matrix – to obtain the pseudo-inverse matrix p and unscramble the FFT. See what comes out as the DC term.
That’s the idea. Could be wrong.
-Bernie
Bernie:
Yours is a fascinating venture into possible full recovery of the bandlimited signal spectrum from non-uniform sampling, which I’ll follow on your site. My objective is far more modest in scope and less demanding analytically: how to reduce the bias of practically estimating the monthly signal mean when only daily Tmax and Tmin are available, as is the case with much of historical data.
The first step is the recognition of the asymmetric wave-form, due to phase-locked harmonics of the typical diurnal cycle, as the source of that bias. Physics of the diurnal cycle provides the first theoretical clue about that asymmetry and the magnitude of eta. Assuming cosine-bell insolation over the 12 daylight hours of the equinoxes, we recognize 1/pi of the peak value as the average level under the diurnal power curve. That energy is redistributed in phase by the capacitive system to produce a temperature response that peaks not at noon, but closer to 3pm. There’s also an attendant reduction in the peak temperature, which doesn’t have time enough to respond fully to peak insolation. Thus the first-cut theoretical expectation for eta is greater than 1/pi, but certainly less than 1/2.
Nick:
If I understand you correctly, you seem to think that there’s much to be gained by sampling once a day at two fixed, but unrelated times and resorting to Fourier analysis to establish the magnitude of the diurnal harmonics. But that puts the Nyquist at 1/2 cycle per day, which aliases all of the spectral lines of the diurnal cycle into zero-frequency. Since there are more than two such lines, you’re left without any means of distinguishing their magnitude.
BTW, time-of-reset (usually misnomered TOBS) doesn’t come close to “dominating” the daily readings of extrema in good station-keeping practice. It’s only the foolish practice of resetting the thermometers at times close to typical extrema and uncircumspect acceptance of values identical to the reset temperature that corrupts the data.
Sky
” It’s only the foolish practice of resetting the thermometers at times close to typical extrema”
I don’t think there is much to be gained by sampling in that way, and in any case no-one is likely ever to do it (though there is quite a lot of Australian data at 9am/3pm). I’m just saying that it (or any other irregular but repeated day-to-day pattern) could, if done, be analysed that way.
I did an analysis of the effect of TOBS on a simulated process using 3 years of Boulder CO data here. The time effect is fairly broad; it isn’t just the obvious extrema times that cause a problem. Of course, when done the operators were not expecting the later use of the data; I believe the NWS encouraged late afternoon reading for the next days newspapers. Nowadays, they work out AWS min/max from midnight to midnight.
“resorting to Fourier analysis to establish the magnitude of the diurnal harmonics. But that puts the Nyquist at 1/2 cycle per day”
The Fourier analysis, based on higher resolution for the diurnal cycle, does establish those magnitudes. My arithmetic just sets out the magnitude of the amount aliased to zero. It would give something for every harmonic of the diurnal, but the trig sums are near zero if the sampling is near periodic. So 7am and 6pm, say, would give answers not so different from 7am/7pm, with small effects from the odd harmonics.
Nick:
My point about the corruption of max/min data series by ill-chosen reset times is entirely missed in your reply:
I’m certainly not advocating “sampling in that way,” which has been much too often the foolish practice for reasons of convenience not only in Australia, but also in the USA. Nor am I suggesting that this accounts for the intrinsic difference between the mid-range value and the true mean. (Although even your Boulder results show a typical sharp reduction in offset when the reset times are far removed from times of typical diurnal extrema.) All I’m referring to is the failure to properly register a particular day’s extremum whenever that happens to be INSIDE the range restriction tacitly set by the temperature at an ill-chosen time of reset.
Regarding DFT determination of the offset from “higher resolution of the diurnal cycle,” your expression, b2*cos(w*t1)+b2*cos(w*t2), reveals a critical dependence of each term upon the actual times tn of the extrema. These much-variable times are never known from historical data and require considerable additional effort to determine accurately from higher resolution modern data. While that effort may be worthwhile for academic purposes, it offers no practical advantage over much-more-direct, time-domain determinations of a simple offset.
Long experience at various project sites scattered throughout the globe shows that properly registered monthly mid-range values are very closely reconciled with true means by the linear combination (1-eta)Tmin + etaTmax. That’s far from the case with Karl’s blanket TOBS adjustments, which jack up modern temperatures inordinately.
P.S. On the other hand, if tn are not related to max/min occurrences, but arbitrarily chosen times, please explain how your expression then relates to the offset.
I think that the min/max would be just fine if weather was static.
The problem lies in the fact that a cold front may move in right after noon sinking the temperature.
In such a case you would have a min during the night and a min in the afternoon with a relatively short high.
How would that affect your average temp?
It also discounts precipitation.
It will be blazing hot out here in the desert, but after a thunderstorm moves through the temperature drops like a rock.
Min/Max may work just fine for certain areas at certain times, but the methodology is not accurate enough for betting the future.
Quote from Jim Gorman:
“To calculate an uncertainty of the mean (increased accuracy of a reading) you must measure the SAME thing with the same instrument multiple times. The assumptions are that the errors will be normally distributed and independent. The sharper the peak of the frequency distribution the more likely you are to get the correct answer. ”
*******
And in the pharmaceutical,from which I come, our GMP processes come into play during our method validation studies. Not only do we constantly test that the methods are repeatable but we also test that it is repeatable across instruments and across technicians. The method, the instrument and the technicians all have to be qualified and each have to be validated. Some of these are tested daily. As well a rigorous and well documented maintenance program.
Now someone, please, tell me that these temperature stations adhere to some form of standard testing and maintenance protocol and standardized manufacturing protocols.
I do have to say that I will still be dubious of a singular derived temperature value, regardless of how mathematically valid it might be, as planet Earth is just too varied and too large of body in terms of geology and climates for single value to have any large planetary practical meaning. The mechanics just do not equate to me.
Thanks to Nick Stokes for the write up. I always enjoy reading his math. I am curious as to which program was used to crunch the numbers: Matlab, Excel, Octave, etc.
“I am curious as to which program”
It’s all done in R.
Nick,
I see how your proposed adjustment offsets the estimated mean to more closely match the true mean. I notice that in what’s left there are still significant differences between the May averages in different years. For example, 2010 versus 2017 shows about a 0.26C difference. If we’re trying to resolve temperature trends on the order of tenth’s of a degC per decade or century is this an acceptable variation?
” If we’re trying to resolve temperature trends on the order of tenth’s of a degC “
This is a single site. Those trends are for the global average.
But yes, I calculated the effect of the common component of diurnal variation. That can’t help with inter-year differences.
Hey Nick,
Thanks for your effort and post.
But the situation can be repaired without resort to high frequency sampling.
Looks like in such discussions the following ‘repairing’ pattern emerges:
Phase one: Contrary to denialist attacks there is absolutely no problem with sampling, aliasing or Nyquist. All is perfectly fine!
Phase two: Well, there may be a slight problem with underlying data but do not worry: all those errors will eventually cancel out!
Phase three: Well, there is bit deeper problem with quality of data, indeed due to aliasing, but we have fixes for that! Retrospectively adjusting and infilling data will make all things right! Climate Adjustment Bureau can fix anything for you 😉
It appears to me that we have reached phase number three.
Every now and then, in climate blogging, one hears a refrain that the traditional min/max daily temperature can’t be used because it “violates Nyquist”. In particular, an engineer, William Ward, writes occasionally of this at WUWT
But that’s not what Mr Ward said. He said that averaging daily min/max creates additional errors that need to be included in the known list of errors and uncertainties inherited from the historical temperature records. That’s far cry from the statement that daily min/max cannot be used for any purposes.
May I ask couple of questions with respect to your text?
1. In the table, column one, how actually sampling is done? Say, sampling every hour – does it mean single sample every hour on the hour, or random single sample from within each hour or true mean of all 12 samples per each hour?
2. My understanding of the anomaly is that it is a difference between a constant baseline value (usually calculated from the longer period as 30 years or so) and a current value. So what do you mean then by expected value? Is it value you somehow expect or is it a current value? For densely sampled temperature data we’ve only 10 or 12 years of data so baseline for anomaly would be period of 3 or 5 years, I suppose.
3. For me the main message of your text is that, knowing that there are diurnal oscillations, we can reduce error due to averaging of daily min and max. Does it mean that if I synthesise a subhourly temperature signal that resembles actual daily temperature signal spanning 50 or 100 years, then extract daily min/max (without timestamps associated with them) and send it to you, you will be able recover very good approximation of the original signal and therefore reduce errors due to min/max approach? That would be an interesting exercise.
Note that this is true for sampling at prescribed times. Min/max sampling is something else.
But monthly averages do not use sampling at prescribed times but daily min/max averaged over each month. So, assuming that this ‘quick fix’ is applicable to reduce error due to sparse sampling it looks like monthly averages still require an additional fix.
Paramenter,
“In the table, column one, how actually sampling is done?”
The raw data has 288 readings/day, every 5 minutes. I arrange it so each day starts at midnight. So when I talk of hourly sampling, I list the 12.00, 1am, 2am etc samples.
“So what do you mean then by expected value?”
I’m speaking rather generally there (E+A). It’s just whatever you’d expect based on some known pattern. For global average, it is the average over a base period for that place and month. Here, because of the finer resolution, it is the expected value for the time of day, for that month and place. The simplest way to estimate is just to average over a period. For example, to get the expected value for 3.05pm June in Redding, I could average the 150 daily 3.05 readings for June days of 2010-2014. The idea is just to remove as much as possible of the daily cycle.
But since I didn’t want to make use of 5 minute data, I would take those 150 days but just the hourly data, average the 150 for each hour, and then do an FFT for those 24 hourly averages. That would, if I chose, give me (after inversion) an estimate of the diurnal to any time resolution I wanted.
In fact all I want to do is to subtract that FFT-estimated diurnal from the hourly data, and month-average that. I know that averaging the diurnal harmonics would yield zero.
“Does it mean that if I synthesise a subhourly temperature signal that resembles actual daily temperature signal spanning 50 or 100 years, then extract daily min/max”
No. This analysis is restricted to periodic sampling. You could send me 6am/6pm samples and a few years (say 5) subset of hourly data, and it should recover an improved estimate of the averages for the 5 min set.
Hi Nick,
Thanks for the informative exchange with Paramenter. I know you will be responding to my other reply, so this is not an attempt to rush that. Take your time. I would like to get some thoughts from you about the trend and offset errors presented. What can we learn from that information? The exchange here with Paramenter presents opportunity to further engage, and if some of it overlaps with my other reply, that is ok.
Do I understand correctly that this correction technique you are presenting would not be of any help to improve the “max/min method”? Said another way, are you saying we can’t correct the aliasing we see when we use the max/min method to process the instrumental record for means and trends?
Also, you mention: “In fact all I want to do is to subtract that FFT-estimated diurnal from the hourly data, and month-average that.” Maybe I’m misunderstanding, but why is the algorithm to subtract? Are you using phase information as well as magnitude? While it is generally stated that aliased energy “adds” to the spectrum of the original, the net effect could be subtractive if the phase relationships are such, correct? So, it would seem that the correction would need to know the phase relationships in order to know which sign to use with the correction.
What do you think?
William,
“Do I understand correctly that this correction technique you are presenting would not be of any help to improve the “max/min method”? Said another way, are you saying we can’t correct the aliasing we see when we use the max/min method to process the instrumental record for means and trends?”
The first is true. For various reasons the Nyquist theory simply can’t be applied to samples which aren’t equally spaced, and where you don’t even know the sample time. And so for the second, the discrepancy with max/min can’t be treated as aliasing. But it is still useful to know how well or badly the month could be estimated if periodic sampling were used.
As to the FFT, yes, phase and magnitude are required. But I’ll set out in response to your longer query how FFT isn’t even needed.
Nick,
Thanks for your reply. We are in-sync regarding FFT and inability to improve/correct max/min method.
You commented: “For various reasons the Nyquist theory simply can’t be applied to samples which aren’t equally spaced, and where you don’t even know the sample time. And so for the second, the discrepancy with max/min can’t be treated as aliasing.”
My response: In the history of sampling, there has never been one instance where samples were spaced equally, as there is always imperfection in the real world. This timing imperfection results in measurement error in the sample. If we sample early or late compare to the perfect interval, we measure a value that is correct for that actual instant, but not correct for where the sample needs to fire for perfect reconstruction. So the question for you Nick, is what is the limit allowed on this imperfection before Nyquist “can’t be applied”? By logical deduction, if no limit can be made and supported by theory and mathematics, then Nyquist applies. It is just that sample amplitude error increases proportionally to the timing error – or there is not enough information to even determine the extent of error.
Furthermore, can you agree that 2-samples/day, with sample timing error or with timing lost (or never recorded) fails to meet Nyquist? (We can’t reconstruct with this data). If it fails to meet Nyquist, how is this different from violating Nyquist?
Can you respond head-on to this line of reasoning?
I believe you said along the way that averages and trends are what we are after and the max/min method allow us to get those. Did you have a look at the trend error I showed from using max/min method? What do you say about this?
William,
“So the question for you Nick, is what is the limit allowed on this imperfection before Nyquist “can’t be applied”?”
I mentioned in the last thread the analogy with heterodyning†. With periodic sampling, that is locked to the diurnal frequencies, producing residues with zero frequency. With near periodic frequencies, you can indeed get aliasing to low frequencies; it’s like AM sidebands. Very low fr aliases will have a similar effect to zero on the monthly mean. Higher frequencies will have less effect.
But as I said, the main issue is not the jitter, but the fact that the samples are determined by value. One reason why sparse samples have trouble is that the result is dependent on the phase at which they operate. Min/max is, in a way, locked to phase.
“fails to meet Nyquist? (We can’t reconstruct with this data).”
Well, there you go again. No-one is trying to reconstruct. They want to get a monthly average. And I showed how sparse periodic sampling affected that. It’s a limited effect that can be substantially remediated. That is what happens when you “violate Nyquist”. Or when you have limited resolution.
“Did you have a look at the trend error”
Yes. I commented here. There is no evidence of bias. You found for 26 stations, min/max gave a higher trend. I found for 109 stations with >10 yrs data, min/max gave a slightly lower trend. But I do not believe any of this is statistically significant.
† To follow up on the heterodyning idea. As 1sky1 said, sampling is like multiplying by a Dirac comb. That is, by a series of narrow but high pulses. That is a form of mixing, except that the comb has the base frequency and all its harmonics at equal magnitude. Jittering the sampling mixes in another signal with possibly low frequency components.
Hi Nick,
You said: “Well, there you go again. No-one is trying to reconstruct.”
My response: Yep. I have said it many times. Its not about whether you are trying to reconstruct or not – it is about whether you *can* reconstruct or not. If you have not sampled such that you can reconstruct or if you butcher a good set of samples such that you cannot reconstruct then a monthly average calculated with these butchered samples gives you error.
Regarding the trend data, perhaps the word “bias” is loaded in a way you object to? I use bias interchangeably with error. Using the 288-samples trend as reference, we can calculate the error or difference from using max/min generated trend. I have seen max/min give both too much warming and too much cooling, but my 26 randomly selected stations show mostly warming. (1 of the 26 showed cooling). I was not suggesting that the error will always be warming. I have said that the error seen is of similar magnitude to the claimed trends based upon max/min. If we had long term data with higher sample rate perhaps we would show even more warming or maybe even cooling instead of warming. I have said we just don’t know based upon the max/min data. I’m not sure how you determined that the errors are statistically insignificant. It seems climate science has been distilled down to a magic machine that gives you even better results when you increase the amount of erroneous input fed into it.
I don’t see the point of bringing up heterodyning, homodyning, or the dirac comb. This is not needed to explain the effects of clock jitter. No one can provide a limit for jitter because there isn’t one, as long as the samples stay in their period. The locking of phase is also not relevant. What matters is what happens if you reconstruct (or if you can reconstruct). That is what tells you if you are compliant or not.
I don’t think we are going to change each others mind. So we can leave this unresolved in that regard, but I appreciate the discussion. Also, a key point in your post here seemed to infer that the error from max/min method could be corrected. You agree it cannot. If we want a more accurate mean and trend then we need to sample at 24-samples/day to meet most signal requirements or 288-samples/day (or more) to meet engineering requirements to satisfy all possible conditions. While it was neat to see you work out the correction of the low sample rate data, it was sort of a round-trip to nowhere. If you have the higher sample rate data (which is accurate to the analog signal), then why discard much of it (introducing error) to then just go back to the higher rate database to gather an average that you can then use to correct the data back to near (but not quite) where you started?
Hey Nick,
No. This analysis is restricted to periodic sampling.
Thanks for providing more details. What it tells me is something like: if we have access to good set of data we can correct errors in corresponding set of bad data. Well, the trouble is that often for historical records we don’t have access to good data. It’s is not something from nothing, this lets call it ‘sampling adjustment’ cannot help with historical daily min/max records where no such extra information is available. This extra bits of additional information required to build ‘sampling adjustment’ come from firstly need for periodic sampling with known timestamps and secondly from what you call ‘high resolution averaging’. Unfortunately, there are no free lunches.
Nice try anyway. At least you did not try to ignore the problem or hand-waving by invocation that ‘it will all cancels out’. Still, would be nice to see more data from different stations, if that approach gives stable adjustments.
“What it tells me is something like: if we have access to good set of data we can correct errors in corresponding set of bad data. “
Yes. But the key thing is that it isn’t the same set of data.You only need one year, or a few years at most, to get a good estimate of diurnal cycle, which you can then use for all other years. Even a totally uninformed guess, say a sinusoid peaking at 6pm, would probably give an improvement, at least for 1/day sampling. OK, you’d need an estimate of range.
Nick,
In the language of analytical chemistry, “subtracting the blank” has been around since the first early colorimeters, Instruments with liquid in a glass cell, shine light through, pass it through a colour filter, measure how much gets through to give a measure of the colour intensity). The method was troubled by effects like suspensions also absorbing light, so it was not long to the twin cell. The other glass cell held the same contents but without the key colour=making chemica.
Subtracting the blank usually gave a better answer. It could also give a worse one, when the optical cells, were not matched, when the photomultipliers were not matched, when one cell got a dirty fingerprint on it, etc etc. Increase the complexity and you increase the error, sometimes.
Nick, you seem to be saying that conversion from absolutes to anomalies improves accuracy. Like subtracting the blank cell reading makes the twin cell device more accurate. It might, but it cannot be done without assumptions.
If you are going to derive a generic curve of known natural T variation for subtraction from observations, you have to be able to show that it is generic, that is, accurate in all times and places. I do not think we can do that, on present knowledge. Certainly we cannot go back decades in history with the same subtraction, because we do not know if it has changed shape or stayed constant. We have lost the past window of time to re-create such curves and we have to assume too much about uniformitarianism if we proceed.
For example, the curve for subtraction could well have past physical variation depending on the incoming solar spectrum and its various components of UV, Vis, IR. We know these varied in the past. We do not know the detail of variation sufficient to construct a generic curve. Besides, in terms of adjustments/corrections, there are bigger fish to fry than those from sampling frequency.
While I agree with the concepts behind your “subtracting the blank” exercise, it might have application over the past few years of satellite observations of a number of physical properties able to affect temperatures, but I cannot see value in application to historic data. As for future data, we now have better gear to “sample” to a pre-determined standard, so we should not need such corrections. The discussion then becomes more academic than practical, at which time I usually retire from it. Geoff.
Geoff,
“you have to be able to show that it is generic, that is, accurate in all times and places”
I’m sure that this process is used in many situations; I’ve used it myself. It probably has a name (or several), but I thought it was worth setting it out explicitly.
It doesn’t have to be accurate in all times and places; it just has to offer an improvement. In this case, most of the aliased error comes from the interaction of sample frequency and diurnal cycle, since the two are locked in ratio. Any E that subtracts out part of the diurnal cycle will be an improvement, as long as E itself can be averaged in a way that avoids the coarse sampling aliasing.
Hey Nick,
Yes. But the key thing is that it isn’t the same set of data.
For me the key message is that anomalies themselves cannot help unless they are supported by good quality data which serves as reference for any adjustment. There may be ways of assuming sinusoidal cycle and from that trying to reduce error but as you said we’re entering here more into guesswork territory – one adjustment build on underlying adjustments.
Nick,
Where you assert about the curve to be subtracted, it does not have to be accurate, it just has to offer an improvement.
This is where your theoretical path departs from my applied path.
I would say you are wrong, because you have no way to discern an improvement.
The best you can do is compare to expert preconceptions, the final retreat method of the IPCC and its departure from proper science.
Interesting discussion, though. Shame it goes nowhere because the magnitude of all errors combined makes much of the historic data unfit for purpose. (Per BOM, for example, the metrication/rounding error alone might be one eighth of our century long global warming, but we have no way to deal with it. And so on into the night.) Geoff
Stokes
You said, “It doesn’t have to be accurate in all times and places; it just has to offer an improvement.” But, unless you can offer an assurance that it is universally applicable, you run the risk of improving some and making others worse, for no net gain.
Nick Stokes, thank you for the essay.
Mr Layman here.
Past temps have been measured as daily Max and Min because that is what the thermometers of the day recorded. The “average” for the day was them divided by 2?
Of course, dividing by 2 doesn’t really tell you what average temperature for the day was because it gives no clue as to how many hours the temps were near the Max or Min.
(For the moment lets ignore what happened after Hansen got his hands on the records.)
We now have systems in place that can record more than the Max and Min but minute by minute temps.
Has there been any effort to look at those current records to see if there is a better number to use than “2” to use for past average daily temps?
(Still huge error bars and still not “global” by any means.)
“Has there been any effort to look at those current records to see if there is a better number to use than “2” to use for past average daily temps?”
First, it is actually usual to record monthly average of Tmax and Tmin separately, and people often get TAVG just by averaging those monthly averages. And then TAVG is what it is. I don’t think anyone has tried to combine in different proportions. The usual focus is on time of min/max reading, which is the basis of TOBS.
We have spent a large amount of time analysing the sampling from individual stations. The biggest problem is the two dimensional sampling input. An example being the sampling interval of places in the US and Europe compared with the number of stations in the Arctic, Antarctic, Africa, Asia etc etc. When you find met offices using predictive algorithms to provide data input for a sequence, then you are on very dangerous territory. An example being the Arctic where there is interpolation used on 1 sample in a 1000 km radius. With such low quality data use of the term global average is meaningless.
Nick:
Glad to see you have the time to dispel all-too-common, amateurish misconceptions about frequency aliasing. It’s astonishing how many readers fail to realize that discrete sampling in the DSP sense is not the sampling of ordinary statistics. Your label “per hour,” however, is misleading in your tables. It’s the periodic sampling interval that you seem to list
Sky,
Thanks. Yes, the heading should be as you say.
As temperature is being used as a proxy for the energy in the system then (min+max)/2 suffers from more than Nyquist as just a sampling error in magnitude for temperature.
Nyquist DOES apply as in spacial terms in the same way as in the size of the cells in any model effect the results derived. If you sample on an unequal spacing then some underlying assumptions are created as to the actual field being measured. This is especially so because of the time slew that exists as the time of day moves across the weather patterns that are in themselves moving also.
The biggest assumption is the peak to peak temperature measurement being a good proxy for the thermal energy under the curve over the same period. That and the fact that 2m is well within the turbulent vertical boundary layer that exists all over the globe.
Depending of the actual profile of the curve over the day, a widely different figure will obtain for the (min +max)/2 and the ‘true’ energy within the system.
I wanted to post an image, but I think I made a mess of it. Administrators may delete my comment.
What’s the right way to proceed for posting images?
Hey Geoff,
Nick, you seem to be saying that conversion from absolutes to anomalies improves accuracy.
My understand is bit different – anomalies calculated directly from the original data cannot help as they simply inherit errors associated with daily min/max. Nick is able to reduce errors by using additional sources of information, not available for historical records. Then and only then we can make couple of not always warranted assumptions about historical records, as that daily temperature profile known from recent years from high quality record is exactly the same as in past years, and apply adjustments.
Well, I’m just a random idiot, but it would seem to me that the rotation of the planet has been happening at a fairly constant speed, and that this would be the dominant factor in this cycle.
Hey Clyde,
But, unless you can offer an assurance that it is universally applicable, you run the risk of improving some and making others worse, for no net gain.
For me Nick claim is far more modest: he’s saying that having access to high quality data we can downsample it to 2 periodic samples per day and reduce the error knowing spectral characteristic of the high-quality signal. For me – it is all. Nick between the verses may be suggesting that we may think how to apply this adjustment with larger scale but this adjustment does not apply when we have only daily min/max from historical records.