GIStemp F to C convert issues

GIStemp, Precision, and FORTRAN

This posting is a bit esoteric, but this stuff really matters. Computer security and computer forensics depend on this kind of detail code review, as does proving a system is reliable and not making really bad errors. Almost as bad, you can have a tiny error, like 1/10% in a line of code that may not matter, but if you have 1000 such lines in a system, you now have a 100% error (assuming the errors add). So “little things can add up”, and computer professionals know that you must use a very very fine comb in a code review to prove a bit of code is ‘valid’.

In looking at GIStemp, one bit of code has regularly been a bit “problematic”:

That is USHCN2v2.f

(Where there was at one time a “bug” that caused it to miss read the input data and drop the 1/100th degree position). This code turns US data in F into “GHCN” format data in C.

Now this data comes from NOAA already a bit “buggered” in that it is in 1/100 F precision. Since the raw data are in whole degrees F, this can only be the result of adding one month worth of those integers, then dividing by integer days, and assuming that the 1/100 F place means something. It doesn’t. Proper behaviour is to either truncate or round to the nearest integer (NINT in FORTRAN). But GIStemp runs with it in 1/100 F of False Precision anyway. OK, I’ve covered that in depth here:

Mr. McGuire Would Not Approve

So we’re just going to assume that issue exists and “move on” to the next issue. The next issued is “what happens with those temperatures with 1/100F precision when processed?” What do they do to the Global Average Temperature?

Into The Code

In looking at that bit of code (and doing some testing of it) the bug that dropped the 1/100 F place seems to have been patched… but… There was something about the way temperatures were converted from F to C that bothered me.

The code in question is very short:

if(temp.gt.-99.00) itemp(m)=nint( 50.*(temp-32.)/9 )

This basically says that if the “temp” variable is bigger than the “missing data” flag, set the “Integer temperature variable – itemp” for a given month to be the “nearest integer – nint” to the result from using the standard formula for C to F conversion (5/9 x temp in F -32).

The only odd bit is that there is a “factor of ten” multiplied in so that the 1/10 C place can be put in an integer with an implied decimal point. (One wonders why not just make it a high precision floating point number, but asking “why some programmer does something” too much is not good for your sanity…)

This ought to be evaluated as (temp-32.) first, then multiply that by 50. and then divide the result by 9 (which, being an integer, would have to be “promoted” to a floating point number to do that division).

OK, what’s the problem?

Well, FORTRAN handles numbers a bit “differently” from what people do. And computer math has some “odd things” to keep track of.

First up is “digits of precision” in formulas. When we multiply by 50, we slide the “decimal point” over by one position for that factor of ten. You could multiply by 5, then move the decimal point over one space. So if temp were, oh, 43.1F, we would turn that into 11.1 (by subtracting 32.), then multiply by 5 for 55.5, then slide the decimal over one spot to get 555. ( Something similar happens in the computer, but in base 2 or binary numbers.) Now we divide by 9 and get 61.67 which we then turn into the “nearest int” or 62. This represents 6.2C and would be fine, but for the fact that we had to add some .03 C in the rounding process.

Or we could divide 50. by 9. and get 5.56. Then we could just take 11.1 times that and get 61.76. We still round up but this time only adding .024 C. With a large bunch of numbers, the difference in these two methods can change the result; especially when they are averaged together. And the first way is more sensitive to the data (since we multiply the data by 50, then divide by nearly 10, making a shift back and forth dependent on data size) while the second way calculates a constant that is kept intact and does not shift bits one way, then the other, losing a few in the process. When temperatures range over a couple of digits of value (or more, from 1F to 100), this can change the precision in the results quite a bit.

Another one is that “9” with no decimal after it. That means “Use an INTEGER 9, not a floating point 9.00000”. When dividing an integer by an integer, FORTRAN does not round, it truncates. Your low order bits can be very sensitive to this kind of effect. Everything to the right of the decimal point can be “messed up” if you don’t keep track of exactly what happens. Do we raise that 9 to a float before, or after the division?

But there are other rules, called precidence rules, that decide in exactly what order things get changed. Do we subtract the “32.” first? Then is it divide by 9 (INT) or multiply by 50. next? (the decimal point makes it a floating point number, 50.000000 and “temp” is a floating point, in this case a “real*4” with lots of precision). So at some point that 9 will need to be turned into a “real” to do the math with the “temp-32.” part. Then we “cast” it with nint back into an integer. So what happens to all those bits to the right of the decimal point? EXACTLY? And does this change the value in the 1/10 C place when it is finally repackaged back into an INTEGER?

Well, the fact is we don’t know. You can try to figure this out (the FORTRAN definition tells you what ought to happen), but it’s better to just test it. Sometimes the syntax of a language is a bit vague about what part gets done first, and some times one compiler (the code that turns FORTRAN into something the computer can run) does things one way when another compiler writer did it another. And sometimes compilers, too, have bugs in them. So you test.

At any rate, I decided to compile this code with both f77 and g95 (two different FORTRAN compilers) and make sure nothing changed. But something DID change. The code is sensitive to exactly which compiler you use.

The Test

How sensitive?

I instrumented the code and did a few runs. I also added a slightly better way of doing this same code that is far less sensitive to such minutia.

What I did was to take that “50./9” and make it into “50./9.” (all “real” data types) and put it into a constant called “ftoc”. This not only makes it VERY CLEAR to the compiler that it does this bit first, but it also only does it one time (rather than once for each record processed) and thus saves about 1.5 Million divide operations each time the program is run, so it’s a lot more efficient. My code now says:

if(temp.gt.-99.00) ttemp(m)=nint(ftoc*(temp-32.) )

Then I ran the program and compared the two methods. In this first case, I use the f77 compiler:

/usr/bin/f77 -o ../bin/USHCN2v2.exe ./USHCN2v2.f

then ran it:

replacing USHCN station data in v2.mean by USHCN_noFIL data (Tobs+maxmin adj+SHAPadj+noFIL)
  reformat USHCN to v2.mean format
extracting FILIN data
getting inventory data for v2-IDs
 USHCN data end in  2007
 Error rate:   0.0104666771
 Old Warmer:   14312.
 Old Cooler:   2120.
 Count:        1570127.
finding offset caused by adjustments

I added some “instrumentation” code that also looks at how many records change, which way, and what the effect might be.

So we have about 1% of the records “different” between the two methods. Of those, in 14312 cases the result is a warmer record by 1/10 C while in 2120 cases we get colder by 1/10 C. The net result will be about a 1/1000C warming of the total of all records (but individual months can be warmed more, and the effect after “reference stations” spread things around is, er, unknown.)

But wait, there is more…

When compiled with the g95 compiler we get similar, yet different results. The original code is “compiler sensitive” to some small degree. In all three cases, the results are different. Though “my way” is not sensitive to which compiler you use. Not surprising given that “their way” will be more data dependent with bits at the far right of the number being lost if you have a 100F day and kept if it’s only 3F.

/usr/local/bin/g95 -o ../bin/USHCN2v2.exe ./USHCN2v2.f
#echo “USHCN2v2.f is compiler senstive at run time. 1 digit rounding variation

  reformat USHCN to v2.mean format
extracting FILIN data
getting inventory data for v2-IDs
 USHCN data end in  2007
 Error rate:  0.009524695
 Old Warmer:  12801.
 Old Cooler:  2154.
 Count:       1570129.

At least the number of records that vary has gone down by a smidgeon.

The Meaning of This Is? About 1/2C Possible from Math

So what does this all mean?

It means that the “false precision” of the 1/10 F place and 1/100 place in the raw data from NOAA influence the conversion to C with an impact that will represent about 1/1000C of warming of the data in total from the process alone. Not a lot, but every little bit adds up. This is just One Line of the code…

And further, because it isn’t “determinative”, you have a compiler sentivitiy issue as well. Luckily, in this case the new compiler makes the planet cooler. Who knew? Maybe we could just buy NASA a new compiler and cool the planet that way…

Just to put a finer point on this, I used 5 variations on simple math to calculate the Global Average Temperature from this dataset. (Not the particular GIStemp code above, just different reasonable choices about calculating GAT from the raw GHCN dataset.)

 JAN  FEB  MAR   APR   MAY   JUN   JUL   AUG   SEP   OCT  NOV  DEC ALLDATA
26    38    73   117   156   187   205   200  172   130    79   39  118
2.60 3.80 7.30 11.70 15.60 18.70 20.50 20.00 17.20 13.00 7.90 3.90 11.850
2.00 3.00 7.00 11.00 15.00 18.00 20.00 20.00 17.00 13.00 7.00 3.00 11.333
2.59 3.81 7.14 11.41 15.22 18.31 20.13 19.65 16.83 12.64 7.71 3.91 11.612
2.63 3.90 7.34 11.78 15.64 18.75 20.58 20.09 17.27 13.03 7.97 4.00 11.913
2.63 3.90 7.34 11.78 15.65 18.75 20.57 20.09 17.27 13.03 7.97 4.00 11.913

The “bottom line” is that I can move the global average temperature by about 1/2 a degree C by my choice of when to use integers, floats, reals, and what order to do the calculation. The first line: uses integers all the way with the number representing 1/10 C increments. It is the most accurate. The second line uses “floats” but in a way that is “done carefully” and does a good job of reflecting the first line values.

The third line has a “divide by 10” (notice, no decimal point…) so it ends up truncating the temperature series in the fractional part. This is clearly losing the added precision of the data (in the full degree F actual records) that could be partialy preserved with 1/10C records.

The fourth line has a bit less impact from a “divide by 10” rather than “divide by 10.” in a different part of the processing. (Adjust each temperature as read, or adjust them at the end in the average. WHEN do you slide the bits over and lose a few?..)

Finally we have the 5th and 6th lines. These create more warming by the simple expedient of doing more floating point math and less integer math (thus preserving more of the false precision through the calculations and letting them have more impact.)

How I’ll Measure GIStemp

For my future tests of the way GIStemp changes the temperature, I’ll be using the method from line 2 as my benchmark. It is, IMHO, the most honest. What does GIStemp itself do? Well, that’s what this whole thing is all about, now isn’t it?

But the bottom line is very clear:

Before we get all excited about changes of 1/10 C (and before we even think about anything in the 1/100 C place) we need to know how much of that change is real, and how much comes from decisions made in the programming…

To give you an idea what the actual difference in the code looks like here are some quoted snippets of my “testing” code:

          icount(m)=icount(m)+1
          count(m)=count(m)+1.

          itmptot(m)=itmptot(m)+itmp(m)
          ttmp(m)=ttmp(m)+(itmp(m)/10)
          wtmp(m)=wtmp(m)+itmp(m)
          rtmp(m)=rtmp(m)+(itmp(m)/10.)

So here I can keep a counter as in integer (icount) or as a float (count).

Then when calculating a running total of accumulated degrees, I can do it as integers in 10ths of a degree (itmptot), as floating point (wtmp), I can turn that itmp from 1/10ths C into whole degrees with a divide (using 10 or 10. causing different results – ttmp and rtmp).

Later in the code, where I create the final monthly averages, I can do it in a couple of ways too:

        itmpavg(m)=(itmptot(m)/icount(m))
        tmpavg(m)=(itmptot(m)/icount(m))/10.
        trnkavg(m)=(itmptot(m)/icount(m))/10
        ttmp(m)=ttmp(m)/count(m)
        wtmp(m)=(wtmp(m)/count(m))/10.
        rtmp(m)=rtmp(m)/count(m)

        igav=igav+itmpavg(m)
        gavg=gavg+tmpavg(m)
        trnkgavg=trnkgavg+trnkavg(m
        tgav=tgav+ttmp(m)
        wgav=wgav+wtmp(m)
        rgav=rgav+rtmp(m)

These are just 6 of the ways of figuring the Monthly Global Average Temperature and then 6 ways of figuring the TOTAL of all time Global Average Temperature. If we did not divide the individual temperatures by 10 (or 10.) we do it here with the total. And the resultant averages are different. Notice that we could also have chosen to divide the total by 10 THEN divide by the count. That, too, would change things in some cases.

So which do you choose? Convert early or late? In floating point math or integer? Why? Which is “right”?

If you think these things all look maddeningly similar, you are right.

Can you really spot the impact of

wtmp(m)=(wtmp(m)/icount(m))/10
vs
wtmp(m)=(wtmp(m)/count(m))/10.

And yet that is the stuff that a 0.64C warmer January is made of… These detailed coding decisions by some random programmer making what they think are reasonable choices can (and do, and in fact must: you can not avoid making SOME decision in how to write the code…) change the results in the 1/10 C place.

And that is why I refuse to get excited about a computer model showing any changes in the temperature data to the right of the decimal point…

In Conclusion

So here we have a specific line of the GIStemp code that warms about 1% of the data records by 1/10C in the conversion from F to C.

We also have a NOAA supplied data set with fractional degrees of F when we know that the fractional part is a lie (or, to be charitable, an error… of False Precision). NOAA calculate a monthly average by a method that may, or may not, be valid and give us data that are not properly vetted for precision. It has a 1/100F precision when no such precision is possible.

And I've done a test on the GHCN data (already in degrees C so we're not going through the F to C conversion "issue") to see just how much sensitivity there is in the "Global Average Temperature" to the kinds of decisions that programmers had to make at NOAA and GIS in writing their code. What I found was that about 0.5 C is “easy” to change, and that about 0.1 C TO 0.2 C is very hard to avoid. So it is almost certain that the 1/10C place owes it’s values as much to programmers and compilers as it does to any actual physical changes in the world.

Any statement about “Global Average Temperature” changing in the tenths or hundredths place is “dancing in the error band” of the calculations. There is no “truth” in any digits other than the full degrees. The whole body of code must be “vetted” and tested in a manner similar to what I did above before there can be any trust given to the “bits to the right of the decimal point”.

Oh, and I’m changing my copy of GIStemp to use the “ftoc” constant method of calculating F to C conversions. It is, IMHO, much better. In doing tests of trends in the data, I will be using the #2 method for my baseline (“careful floating point – true to the integers”) where possible. It, too, is better.