Wednesday, July 17, 2019

Comparison between global temperature indices following GHCN V4; changes since 2015

I had noticed that recently the concordance of GISS with the more advanced TempLS methods seemed to have improved, and I wondered whether there might be a general improvement associated with the adoption of GHCN V4, with the big increase in land stations. In 2015, I posted a study of the extent to which a rather large set of indices mutually agreed. It included land, SST and troposphere measures. I may revisit that. But for the moment, I want to look in a similar way at just the surface (land and ocean) indices. Since they seek to measure the same thing, differences can be attributed to method rather than physics.

In that earlier post, my measure was the standard deviation (sd) of differences between monthly index values over the most recent 35 years. That was to fit with the satellite data, which is not used here. But I will stick with the period (updated), because while there doesn't seem to be much sensitivity to the choice, I want to concentrate on method differences rather than data, which might diverge at longer times ago. Data sources are listed here. The sd measure is not affected by different anomaly bases.

So here, as an overview, is the current set of standard deviations, according to the color scale in the key on the right. Red means better agreement.

The best agreement is between the various methods of TempLS, as described most recently here, with an overview here. It is so much better that I have used a color out of the rainbow scale to show it. Differences of the three advanced methods have sd's of about 0.01°C. That is about a third of the nearest difference between indices from different sources.

The next best agreement (0.027°C) is between TempLS grid and NOAA Land/Ocean. I commented,also in 2015, on how NOAA and TempLS grid were eerily close; I showed comparison graphs. That closeness has persisted, and is a reason why I keep posting TempLS grid, which I otherwise think is a very primitive method. So the fact that NOAA is so close makes me worry about that index. But anyway, it is now even closer. I have modified the grid method to use a cubed sphere mesh, which I think is much better than lat/lon.

Almost as good is the agreement between GISS and the advanced TempLS methods. As I shall show, this has improved with V4. TempLS LOESS and Infill have lowest sd, at about 0.031; mesh is a little more at 0.035°C.

The five non-TempLS indices are shown in the top corner. Their levels of agreement are much lower. The Cowtan and Way kriging index has an sd of 0.45 with both GISS and BEST, but less agreement with NOAA and HADCRUT. The best agreement (0.039) is between HADCRUT and NOAA; these have always seemed to act as a separate grouping. GISS and BEST agree about as well (0.045) as the do with C&W. BEST has the greatest disagreement, with both NOAA and HADCRUT.

I posted the data back in 2015, so now I'll use it to show how these concordances have changed. In the following plot, the current sd is divided by the sd reported in 2015. A red value indicates reducing sd (improvement). TempLS LOESS is omitted because it did not exist in 2015.

The biggest changes are associated with TempLS, where methods have improved, particularly with Infill. In 2015 this was a heuristic method, which seemed to give a large improvement. But now I solve a diffusion equation to convergence, which seems to be better again. The sd with GISS is about halved, and is, by a hair, the best agreeing of any TempLS. Because it shifts further towards the other advanced TempLS methods, it moves away from the grid method, and so also from NOAA, which shows as decreasing agreement. The improved agreement (4x) with TempLS mesh is the greatest change of all.

The other marked changes are with BEST. 2015 was still fairly early in its life cycle, and most noticeable is the increasing disagreement with NOAA and HADCRUT. But it also doesn't agree with anything very well.

The other indices, interacting with each other and with TempLS mesh, show little change. T mesh was stable over that period. There is some deterioration of agreement between HADCRUT and GISS, which could be due to the introduction of ERSST 4 and 5, which adjust for the introduction of drifter buoys in SST measurement. HADSST is just bringing out V4 which may implement that.

Here is a more detailed quantification of the changes. There are 9 plots, showing for each index the sd's of the differences with the others (green). Overlaid in transparent blue is the corresponding sd from 2015. For TempLS LOESS, I have used the 2015 sd's of TempLS mesh. Use the arrows below to cycle through the plots.
In the first plot (GISSlo) the TempLS advanced indices (TM, TL, TI) show best agreement, and also improvement (faint blue is 2015). Agreement with HADCRUT is worse. Of the other plots:
  • HA HADCRUT - almost everything is worse, especially BEST. The best agreement is with NOAA and TempLS grid.
  • NO NOAA - not much change, except for lower agreement with BEST. But not a high level of agreement.
  • BE BESTlo - again much increased, and high, disagreement with NOAA/HADCRUT. Otherwise small changes toward more agreement.
  • CW Cowtan and Way - much improved agreement with TempLS; fair agreement unchanged elsewhere.
  • TM TempLS mesh - good and improved agreement with GISS and TempLS grid. Very good agreement with TempLS LOESS and Infill, with Infill much improved (due to Infill method improvement).
  • TL TempLS LOESS - as for mesh. LOESS did not exist in 2015.
  • TI TempLS Infill - very good and improved agreement with TM and TL. Also improved wrt GISS and CW; HA, NO, TG somewhat worse.
  • TG TempLS grid - mostly substantially improved, and not bad, except for BEST and CW. Slightly worse relative to BEST and TI. The good, and further improved, agreement with NOAA has been noted.
Overall, I think it is important to note that even the worse disagreements are not so bad - about 0.075°C. There is a marked tendency to clump, with HADCRUT/NOAA/TempLS grid as one group, and GISS+TempLS(TM, TI, TL) as another, with BEST and CW more loosely attached.


To put the size of these differences in context, they range from 0.01, which I called very good, to about 0.075, which was about the worst. But I did a quick similar analysis between HADCRUT, UAH and RSS. The result is here:

The best agreement there is between the satellite measures, as about 0.1°C. Agreement between surface and satellite is in the range 0.125 to 0.145°C


I have posted the data for this post on a zipfile, with readme.txt, here.

Tuesday, July 16, 2019

GISS June global up 0.07°C from May.

The GISS V4 land/ocean temperature anomaly rose 0.07°C in June. The anomaly average was 0.93°C, up from May 0.86°C. It compared with a 0.067°C rise in TempLS V4 mesh. As with TempLS, it was the warmest June in the record, also by a considerable margin (0.11°C).

The overall pattern was similar to that in TempLS. Hot in most of Europe, extending into Africa and the Middle East. Cool in W Siberia, but hot in the NE. Cool in US/Canada, but warm in S America. Antarctica was mostly cool.

As usual here, I will compare the GISS and earlier TempLS plots below the jump.

Wednesday, July 10, 2019

Revised June global surface TempLS up 0.067°C from May.

I made an error in the previously posted TempLS for June. The rise is now 0.067°C instead of 0,096°C. The June anomaly was 0.782°C instead of 0.811°C. The difference in global average is actually small, and 2019 was still easily the warmest June in the record, but by 0.07°C, not 0.1°C.

Although the overall difference is small, the error was major - my calculation used May SST values, not June. Teething problems with the V4 system. Fortunately, the hemisphere effects more or less balanced, but the map looks different. I had noted  earlier that there was a marked hemisphere difference. I did look further into that, which revealed the error, but I should have twigged sooner. Land is unaffected, and so most of my previous comments still hold.

Here is the revised temperature map, using the LOESS-based map of anomalies.



The original map is preserved on my tweet here.

Saturday, July 6, 2019

June global surface TempLS up 0.096°C from May.

Update - see revision here. The average is down by just 0.03°C, but SSTs were wrong month, and the map now looks different (correct version below).

The TempLS mesh anomaly (1961-90 base) was 0.811deg;C in June vs 0.715°C in May. This contrasted with the drop (0.056) in the NCEP/NCAR reanalysis base index. It was the hottest June in the record, 0.1°C higher than June 2016.

I am now showing TempLS LOESS as the alternative (rather than grid); I think it is about as good a method as mesh. It showed a rise of 0.114°C.

There was a marked global pattern (caused by error - see update), with tropics and SH mostly warm, and the extratropical NH cool, with the notable exception of Europe, which was very warm indeed, and NE Siberia likewise. The mostly cool Antarctica was also an exception.
Here is the temperature map, using the LOESS-based map of anomalies.



And here is the map of stations reporting:




Wednesday, July 3, 2019

June NCEP/NCAR global surface anomaly down 0.056°C from May

The Moyhu NCEP/NCAR index fell from 0.41°C in May to 0.354°C in June, on a 1994-2013 anomaly base. That is the third successive fall since the high point in March, and brings the temperature back to between January and February. Still warmer than June 2018 or 2017, and close to June 2016.

It was mostly cool in N America, warm in Europe (except near Atlantic). Quite warm in NE Siberia/Alaska, and mixed, but mostly cold, in Antarctica.

The BoM ENSO Outlook has been reset from Watch:
"The ENSO Outlook has been reset to INACTIVE. The immediate likelihood of El Nino developing has passed with ENSO-neutral the most likely scenario through the southern winter and spring."




Tuesday, July 2, 2019

Fake charge of "tampering" in GISS

I was commenting on an interesting post (part of a series) at Clive Best's blog. He's been looking at the differences between Hadcrut 3, of about 2012 vintage, and current Hadcrut 4.6. There are some, and I may blog about that. The most obvious difference is that the number of stations in the inventory has nearly doubled. But Clive was focussing on changes to locations that were common to both. I did some analysis, part reported here.

As is apt to happen, there were undercurrents that data is being manipulated for some underhand purpose, and Clive was entertaining the idea that the Pause was being suppressed. Not jumping to conclusions, though, but some were more inclined to. There has indeed been a noticeable increase over those years in the trend during the Pause period. This is overdue, since Cowtan and Way showed in 2013 that HADCRUT's deficiency in Arctic stations was responsible for the difference in Pause trend between theirs and other indices.

Anyway, among dark talk about Hadcrut suppressing the Pause, Paul Matthews commented that GISS had done the same thing, and between 2017 and 2019. This surprised me, because I follow GISS, and compare it with TempLS, and did not know of such changes, which if present would presumably relate to transition from GHCN V3 to V4. Gavin Schmidt has also said that the effect of this was very small.

So I followed Paul's link, which led to a Tony Heller post titled "Tampering Past The Tipping Point". It showed the following plot (followed by many more):



And as usual there, the plot and post seem to have circulated widely. You can see a long Twitter listing here of tweets linking to it. So what is it based on?

As often with Heller's posts, it isn't about what most of his audience thinks it is, but they don't seem to worry about fine points. It isn't the GISS land/ocean (LOTI) that gets widely circulated and discussed. The heading says "GISS Global Land Surface anomaly". But GISS doesn't have a Land Surface anomaly index, unlike NOAA or HADCRUT (CRUTEM). So my first thought was that he was plotting the "Met Stations Only" index, Ts. He has done that before, and the years quoted (2000 and 2017) do correspond, more or less, to what is supplied on the GISS History Page (scroll down to where "Met Stations" appears in the headings). I'll digress a little to explain this index.

GISS Ts index

GISS Ts is no longer shown on the main page, although it did have more prominence in V3. Now it is relegated to the History Page, with the introduction:
"For historical reasons we also maintain a calculation of the anomalies that would result if one only used the meteorological station data. This estimate is not affected by issues in ocean data processing, but because the land is warming faster than the ocean, it has a larger trend than the land-ocean index that is now our standard product. That too has been remarkably stable over the years:"
And with that, they give, as they do with LOTI, a plot of the data as it had been presented at various stages of GISS history, going back in fact to 1981. You can see both plots of the curves together, and their differences from current. And indeed the differences are small, especially recently.

The "historical reasons" are that, until about 1995, there didn't exist a dataset of sea temperatures of anything like the duration of the land record. So when Hansen and Lebedeff in 1987 published the ancestor of the GISS index, they used whatever station data they could get to estimate surface temperature over the oceans as well as land. Islands had a big role there. This index, called Ts, or GLB.Ts, was their main product until the mid '90's, when it was gradually supplanted by LOTI, using ocean sea surface temperatures (SST) as needed, as they became available backward in time.

Update. As CCE notes in comments, with GISS V4, the Ts index is not only relegated to the History page; it is not calculated in V4 at all. The numbers I have used are the latest V3.

GISS Land

However, Paul insisted that there was a land index, and pointed to the Analysis Graphs and Plots page. If you scroll down to the heading "Annual Mean Temperature Change over Land and over Ocean" and open, it shows a plot of anomalies over land and over ocean, and below it gives links to data.

Now this is something different to GISS Ts. It also uses station data, but to estimate the average for land only. All such averages are area-weighted, but here is is just by land area. So from being very heavily weighted, island stations virtually disappear, since they represent little land. And the weighting of coastal stations is much diminished, since they too in Ts were weighted to represent big areas of sea.

The important message here is that Ts and Land are not the same, which I will now show with some graphs. Data is sourced and linked at the bottom.

Recent History, 2017 and 2019

Tony Heller provided a spreadsheet with his post, and it had the GISS data for versions of Ts up to 2017, and the Land data for 2019. I have described details of this here and following. But GISS Ts does of course go to present (May 2019), which is regularly posted here. And you can get past versions of the Land average plots with data on the Wayback Machine - here is version of Jan 2017. So let's look at annual Ts, with 5 year running smoothing:



They are actually very similar. I'll givea combined difference plot later. What about Land?



Not quite as close, but also similar. The main difference is that pre-1900 is warmer in the current version, reducing the trend since 1880 from 1.05 °C/century to 1.0 °C/century. The trend of Ts also reduced slightly. Not much sign of data tampering here! In fact, given the number of extra stations in GHCN V4, there is remarkably little change.

Now I'll plot the Ts and Land averages superimposed on Tony Heller's "tampering" plot. But because the 2017 and 2019 versions are so similar, the plot gets cluttered. To make better use of space, I have truncated some of the big colorful annotations. I'll plot just the 2017 version of Ts and the 2019 version of Land. Not coincidentally, these are the versions of each found in Tony's spreadsheet.



They superimpose exactly! What has been presented as a "tampering" is in fact a plot of two different datasets, representing two different things. To emphasise that, I'll now plot 2019 versions of both Land and Ts:



Also a very good fit. The difference between the red and the green curve isn't "tampering" over time. It's the same difference if you take the current versions. They are just two different datasets representing two different things.

Getting it right.

As mentioned, I originally set this out in comments at Clive Best's site, where Paul Matthews first raised the Tony Heller post. I then noted that at that (Heller's) site, a commenter Genava had observed that the 2019 data plotted was different from the 2019 Ts data, which was the index of the 2001 and 2017 versions. That was on June 27. It got no response until Paul, probably prompted by my mention, said that the 2019 data was current Land data. I don't think he appreciated the difference between Land and Ts, so I commented June 28 to try to explain, as above. Apart from a bit of routine abuse, that is where it stands. No-one seems to want to figure out what is really plotted, and comments have dried up. Meanwhile the Twitter thread castigating "tampering" just continues to grow.

Data

The data plotted are year versions of the GISS Ts Met Stations Only index and the GISS annual data for Land Averages. The sources are, in ascii format:
GISS T2 current (2019) version
GISS T2 historic, includes 2017 version in zip file
Land average current, csv format
Land Average 2017 wayback version, txt

The data I used are in a .csv file here.

Tuesday, June 18, 2019

GISS May global down 0.11°C from April.

The GISS land/ocean temperature anomaly fell 0.12°C in May. The anomaly average was 0.86°C, down from April 0.97°C. It compared with a 0.149°C fall in TempLS V4 mesh. I should note that there were late moves in both indices. GISS' April figure dropped 0.02°C between posting last month and now, making a total drop of 0.13. I note that GISS says that this month is the first using GHCN V4 (other than beta); the change may be part of that. The drop in TempLS is somewhat less than originally posted. I posted that rather early, since all major countries seemed to have coverage, and I got stable results two days in a row. The latter proved illusory, though. A minor nuisance with GHCN V4 is that they keep adding to the inventory. I have resolved to only update it once a month, and this was the first time it had been updated automatically. It didn't quite work; the net result was that my apparent stability was due to data not being updated. It should be OK in future.

The overall pattern was similar to that in TempLS. Warm in most of Asia and Africa. Cool in central Europe and much of the N America (but warm in NW). Arctic warm.

As usual here, I will compare the GISS and previous TempLS plots below the jump. I have now switched from my spherical harmonics based TempLS plots to LOESS based, which gives a more similar resolution.

Thursday, June 6, 2019

May global surface TempLS down 0.177°C from April.

The TempLS mesh anomaly (1961-90 base) was 0.703deg;C in May vs 0.88°C in April. This was a lot more than the drop (0.086) in the NCEP/NCAR reanalysis base index, and takes it back to January's temperature.

The drop was mainly in land temperatures, and central N America was quite cold, as was central Europe, and less so, China/Mongolia. Africa and Brazil were warm, and also the Arctic. SST has been drifting down slowly, but is still warm.

Here is the temperature map, using now a LOESS-based map of anomalies.



And here is the map of stations reporting:




Monday, June 3, 2019

May NCEP/NCAR global surface anomaly down 0.086°C from April

The Moyhu NCEP/NCAR index fell from 0.496°C in April to 0.41°C in May, on a 1994-2013 anomaly base. It was a substantial fall, exactly the same as from March to April, and yet it still leaves May as warmer than February, or any earlier month back to May 2017.

The Arctic was quite warm, including Greenland and west. There was a band of cool from E Canada to California. Also W Europe was quite cool, and also China through to central Asia. Antarctica was mixed.

The BoM ENSO Outlook is back to Watch - still an enhanced likelihood.





Friday, May 31, 2019

GHCN V4 Monthly temperature data displayed on an active sphere

This is another in my series on a close examination of the GHCN V4 database of station monthly average temperatures. The previous post described a system for displaying time series graphs and numerical data for named stations (with search). This post draws attention to an older system, where a triangle mesh is shown for each month, with the reporting stations as nodes, and shaded according to temperature anomaly (TAVG) for that month. The shading is such that it is correct for each node (with complications where some are nearly coincident), so it reveals a lot about the low-level structure of he data. There is enough homogeneity that the overall pattern shows through. I had thought about replacing it with a LOESS-based map, which would actually be a lot easier, and is better for the overall picture. I'll probably put that up at some stage, but I think the low-level info of this plot justifies its place.

The page for this facility (link on right - "WebGL Map of...") is ancient, more or less for the duration of GHCN V3. The original post from 2012 is here. The plot is now based on the Moyhu WebGL facility. That means that as well as being able to drag the plot around as a trackball, you can click on the small map for quick centering. The overlay of mesh and points can be toggled off and on with the checkboxes beside the small map. The table beside that lets you choose a month. With the top bar you choose the decade, next the year, and then the month. When the right date shows below, click "Fetch and Show". When the right data is showing, the date will be on top left. There are now zoom buttons, as well as zooming but dragging vertically on the globe with right button down. The Subset button lets you speed up response by choosing a subset, but I don't think you'll need that here. You may notice that I have shifted from Rainbow to GISS-like colors and ranges.

When you click on the sphere, it shows the code name, name and anomaly of the nearest station. The anomalies are as calculated by TempLS V4 mesh, and the base period is 1961-90. The color scale is centered to the mean for the month, so it varies as you go back in time.

I'm using unadjusted data, and there are a number of stations that stick out as different. I actually developed the viewer to see what is going on there, and I'll write about that next.





Wednesday, May 29, 2019

Viewer for GHCN V4, with search and active graph

I have made an active viewer for the GHCN V4 database of station monthly average temperatures, as used in TempLS and other indices. I was prompted to do so after adapting the GHCN webGL page to use GHCN V4. I think the results are generally good, and I'll write more about that soon. But there are a few stations that stand out as having very odd anomalies. The viewer shows adjusted and unadjusted temperature data for each of the 27380 stations in current GHCN V4. I'll link the page more permanently through the data portal page on the right panel.

The main viewing mechanism is an active graph, from which you can download station data. The pink rectangle below the graph (below) has a search box. You can put in a part of the station name. The GHCN 4 convention is that stations are in all caps, with underscores replacing spaces. The search will return up to 20 names; each line giving the full name (to 18 chars), the GHCNv4 station number, and the approximate year range of the data. To the left of each line is a radio button. Clicking this will download the data for that station and show an active graph.

Below the graph is a table of checkboxes with names of months. There is a column of unadjusted, and one of adjusted. Check the boxes to determine which graphs are shown. There is a legend box below the graph, which shows the selected names and color. You can drag this onto the polot to act as a legend.

The graph is active (like this one). That means that you can drag with mouse the location of the plot. If you swipe horizontally below the x-axis it will shrink or expand the x-range, and similarly behind the y-axis. So you can adjust the plot to see best what you want.
Update If you don't like a color you can change it. Just click the button of that color in the legend to randomly choose from a palette of 64.


The Data button above the checkboxes creates a new window with a printout of the data that is showing. You can save this, or copy and paste.

I should note that I have excluded any data with a quality flag; sometimes this is quite a lot. Anyway, here is the plot:

Tuesday, May 21, 2019

Higher resolution graphics for monthly surface anomaly.

Update: Following suggestions in comments, I have made a new tableau in which the new LOESS plots are compared not only with MERRA but with the ECWMF reanalysis ERA5. I think that where LOESS and MERRA disagree, ERA5 is often in the direction of LOESS. But you be the judge. There are the same 12 months of 2016 in the set.

I have for several years posted a monthly plot of global surface temperature anomalies, calculated by TempLS, and followed up a week or so later comparing it with the corresponding plot from GISS NASA. I use the same style. I chose GISS because it is, IMO, the best of the major organisation plots.

The TempLS plot uses spherical harmonics (SH) to give a smoothed stable plot to anomalies defined at scattered points. It has worked well, but the smoothing restricts to contours of a rather limited curvature, and of course restricts resolution. There is also an issue, noted in my previous post, that the spatial distribution of the anomalies calculated on a fixed interval, like GISS 1951-80, does not quite correspond to mine, which are derived on another basis (least squares fitting) and adjusted to the time interval by adding constants that make the means match. Ensuring that correspondence has some minor effects.

I have a new method of averaging using LOESS on a regular array of nodes derived from a regular icosahedron. This also makes possible a plot of the anomalies based on this local averaging which has apparently higher resolution than the SH based plot. The problem with SH is that resolution is limited by the areas of weakest coverage, while LOESS is truly local, and gives good coverage where the coverage is good, which for GHCN V4 is most of the world. I want to show that the apparent higher resolution is real by comparing with NASA'a MERRA reanalysis, which is posted to approx 1/2 ° resolution.

One might ask - why not just use MERRA? I don't present MERRA as the gold standard of accuracy. But in any case, it isn't available with the immediacy of surface data. I got the MERRA data from KNMI Explorer, where it is current up to end 2016. So my comparisons are for that period.

State of the art

The traditional style of mapping anomalies is to just show the lat/lon grid used colored by temperature. Here you see HADCRUT and NOAA:



In fact HADCRUT has been smoothed and is actually not such a bad plot, except for the white blanks. The NOAA plot (actually land only, from here) is mainly cited by contrarians to claim poor coverage. NOAA assists them here by posting a map from early in the month (Jan 11 here) and then not updating as more data comes in. But both have a similar failing. They treat the arbitrary grid segments as telling something about the data. But they would tell something different if a different grid size was chosen. You don't suddenly become ignorant about temperature because an arbitrary grid boundary has been reached. NOAA seems to realise this, because their corresponding land/ocean plot looks like this:



They have done a lot of infilling. I've never seen an explanation of the basis used, but here I mainly want to focus on resolution, which is still poor. The corresponding GISS plots come with two stated interpolation lengths:

On left is 250 km, right 1200km. The left should have higher resolution, but is again marred by grey blocks where it is deemed that information ceases at a cell boundary. The right fills more area with more smoothing, and is the one I use, though it still has some grey.

The new method, compared

Details of the new method are:
  • The anomaly data is first converted to the 1951-80 base, using the distribution described here.
  • An array of 5762 23042 regularly spaced nodes is created on the Earth's surface, based on dividing the regular icosahedron, with a transformation on each triangle to preserve equal area triangles (described here).
  • A first order LOESS (local regression) based on nearest 20 stations is used to interpolate each month's anomaly data onto these nodes (described here)
  • Using the regular icosahedral grid, those data are then interpolated onto the image of a 1° lat/lon array, projected onto the sphere
  • The data thus interpolated is then used to made a 2D lat/lon plot, using, as I usually do, the GISS temperature levels and colors.
I do this for each month of 2016, and compare with
  • The two graphs I usually post - GISS 1200 km and TempLS Spherical Harmonics, both 1951-80 base, with TempLS converted regionally
  • the new method is bottom right, and on right, the MERRA plot (data from KNMI), also converted regionally from 1981-2010 to 1951-80.
With the buttons at the bottom, you can cycle through the months of 2016. December displays first.



The LOESS bottom left plot looks to have more detail than the top plots. See, for example, the US plot. To test whether the appearance of resolution is real, compare with the adjacent MERRA plot. I'll mention these points in the December plot:
  • First the weakest point - LOESS shows a big cold spot in Alaska, which is not really in either MERRA or GISS. However, if you look at the 250km GISS plot above, there is a cold spot there too.
  • For Canada and ConUS, the extra detail of LOESS does seem to align with MERRA
  • The extra detail of the hot spot over NW China seems to line up
  • MERRA has a lot of detail about Antarctica, where GISS is very patchy. TempLS SH is sounder there, but LOESS moves a little more in the direction of MERRA
  • LOESS, like GISS and MERRA, has an El Nino jet, where SH is rather wavery.
You can check for other points here, and in the other months of 2016. But I think in almost all cases, where LOESS has extra detail, it will correspond to something in MERRA.

I plan to use the new graphics scheme for the WebGL page for past monthly anomalies. That page is currently hard to maintain because GHCN V4 keeps introducing new stations, which means the data does not match the stored meshes. LOESS will fix that, as well as being better.

Update - including ERA 5

In comments, Bryan Oz4caster and Zeke Hausfather recommended using the ECWMF reanalysis ERA5, perhaps in preference to MERRA. I certainly think it is a good idea to have another reanalysis, so that where TempLS disagrees with MERRA, which if any is wrong. So I downloaded ERA5 from KNMI, and made another tableau below. This time the two reanalyses are on the right, and I have duplicated TempLS LOESS on the left. For ERA5 I have used 0.5° grid spacing, similar to MERRA.



My general impression is that MERRA sometimes seems overly dramatic, and ERA5 is more in line with LOESS on such occasions. Looking at December, the cool in western Canada is similar between LOESS and ERA5. The blobs of warmth around Mongolia in ERA5 and LOESS are similar; MERRA a little different. But there still isn't any support for the cold spot in Alaska.









Monday, May 20, 2019

Changing anomaly base in spatial plotting.

This post tries to give a more exact treatment to comparing global temperature graphs with different anomaly bases. It is preparatory to one where I try out a new style of graphics which I think has better resolution than the spherical harmonics (SH) based graphs that I use now, and compare with GISS. I want to compare with high resolution graphics from reanalysis, but first there is something I need to improve.

Temperature anomalies are created by subtracting from each reading of temperature an estimate of the normal, or expected value, for that time and location. Various estimates can be used, as long as they are consistent. Often used is a temperature over a three-decade interval. I use a least squares method, but normalise it to zero average over the 1961-90 period. That last step, though, is an adjustment at global level.

Comparing anomalies with different base periods is normally done by aligning the global average over those period. This works well at that level. But it does not align the anomalies in different spatial regions, where the local averages may have changed differently in those periods. In comparing GISS with 1951-80 basis with TempLS, with actually a least squares basis, I have just added an adjustment constant for each month. I set out the process here, giving a table of changes to make for all the popular conversions. But more is needed for spatial alignment, although the actual effect on say a month anomaly plot is small.

Practicalities - spherical harmonics again.

However, it isn't obvious how to do that. The problem with three decade bases is that not all stations have data in those times. That is why I use least squares. To make a comparison between two periods, you would, in principle, have to limit to stations where you can calculate averages in both periods. Even if that were possible, it would be a pity to have to in effect redo two temperature anomaly constructions just to change base.

But the change that matters is what shows on a spatial plot, and rather than do every station, it is sufficient to approximate the difference with spherical harmonics. This is the method I normally use to show TempLS monthly maps. Since the issue is what the regional temperature difference actually is, rather than the individual measurement methods, it is sufficient to work out coefficients by just one method. I naturally use TempLS.

Method and numbers

I first calculate a spherical harmonics fit for the TempLS anomalies for each month since 1900. It is like a Fourier fit, but using regression. The integration method is the mesh method. I specify order 10, which actually gives 121 functions in the SH basis. For each order n, one of the functions is just cos latitude multiplied by sine of n times the longitude, so this gives a measure of the spatial frequency. More details of spherical harmonics are here.

Once the SH approximation for each month is made, there is no further reference to which stations are reporting. Since the operations are linear, for each tri-decade we can just average the 121 coefficients, and then difference to give the change between bases. Then to regenerate the temperature map, it is just necessary to compute the values of the harmonics at whatever points are being used to make the map, for example a lat/lon grid, and then combine these using the difference coefficients.

Results

First, here is a map of the TempLS values for mean April between 1961 and 1990, plotted with negative sign. It would be similar for other months, since it mainly represents changes over a thirty year period.


The mean should be approximately zero, since TempLS was set to have mean zero in this range. It isn't quite zero, because the mean represents the SH-enhanced mesh integral. It may not be obvious from the plot (in GISS intervals and colors) that the mean is zero, but the area of the very warm part is small relative to the large areas of S Hemisphere cool.

This plot gives an idea of the magnitude of differences to expect. They can be large. Now I'll show the practical effect, for the most recent month, April 2019.



The top two maps are those I displayed. Left is the SH representation of normal least squares TempLS, and right is the GISS plot, which is based on actual means of months in 1951-80. So below it I have put the plot of TempLS calculated on this basis, as described above. Bottom left is the map of expected corrections in going from TempLS to a true 1951-80 base, less the spatially constant offset that was used.

There isn't that much difference between top left and bottom right, or GISS for that matter. But the bottom right is closer to GISS in some respects, and that can be explained by the bottom left. The cool spot above NW Canada is much reduced at bottom right, in better agreement with GISS, and that is the effect of the warming correction shown bottom left. The cold of the W Sahara almost disappears, as it does in GISS. Bottom right shows the warm correction there. Same in Saudi Arabia. Finally the hot spot in NE Siberia is made even hotter, in broad agreement with GISS, although the shape there is slightly different. Generally the ocean corrections are too small to notice.

I'll use the corrected versions in new monthly reports. But mainly I am setting this out because I am planning to use what I think will be a better resolution map. I will show that, in the next post, in comparison with reanalysis, and for that I need the extra accuracy.











Friday, May 17, 2019

GISS April global down 0.12°C from March.

The GISS land/ocean temperature anomaly fell 0.12°C in April. The anomaly average was 0.99°C, down from March 1.11°C. It compared with a 0.097°C fall in TempLS V4 mesh. As with TempLS, because March was so warm, dhe drop still meant that last month was the second warmest April in the record.

The overall pattern was similar to that in TempLS. Warm in all of Eastern Asia, especially Siberia. Warm in Europe, NW Canada and Alaska, and much of Africa. Mostly cool in Antarctica.

As usual here, I will compare the GISS and previous TempLS plots below the jump. I am planning to post in the next few days about a better style of TempLS plotting

Saturday, May 11, 2019

April global surface TempLS down 0.097°C from March.

The TempLS mesh anomaly (1961-90 base) was 0.872deg;C in April vs 0.969°C in March. As with the very similar drop (0.086) in the NCEP/NCAR reanalysis base index, that still makes it a very warm April, the second warmest in the record. It was the warmest month since March 2017.

GHCN V4 results were a little later this month. I've noticed a reversion to the V3 custom of reporting by country; earlier it seemed that most countries were represented at an earlier stage. I'm a bit late reporting, too, because I found that the preliminry results changed from day to day. But I think all the main countries are in now.

This time the very warm Arctic region was in NE Siberia, with a rather cool spot NW of Canada. Europe was warm, and China, especially W.  Africa was warm, the Middle East rather cool.
Here is the temperature map, using a spherical harmonics approximation to the mesh results.



And here is the map of stations reporting:




Friday, May 3, 2019

April NCEP/NCAR global surface anomaly down 0.086°C from March

The Moyhu NCEP/NCAR index fell from 0.582°C in March to 0.496°C in April, on a 1994-2013 anomaly base. March was a warm spike, though, so that still leaves April as the second warmest April in the record.

It was very warm in the Arctic, including over NE Siberia and Greenland, but once again cold over the Canadian Archipelago. Cold in Spain and the W Sahara. There was warmth From Mongolia through to N India and Thailand. Quite cool in the Middle East, from Egypt through to Iran.

The BoM ENSO Outlook is remains at Alert - but "...This indicates that if El NiƱo does develop, it is likely to be short-lived and weak." Still, it's already warm. .





Wednesday, April 24, 2019

Tests and results from TempLS v4 global temperature program.

I have been writing a series of articles leading up to the release of V4 of the prgram TempLS, which I run on the monthly data from GHCN (land) and ERSST V5. The stimulus for the new version was the release of GHCN V4. I use the program to compare with the major indices, which are basically GISS, NOAA, HADCRUT and BEST. I'll make comparisons of the TempLS output with those indices here.

I have been developing different integration methods, with the basic idea is that the agreement between a number of good methods with different basis is a guide to the amount of uncertainty that is due to method. I think it is quite small, and I will show comparative graphs. The key thing is that the different methods within TempLS agree with each other better than the indices agree among themselves.

One of the methods I have been exploring is the use of Spherical Harmonics (SH). This is not so much an integration method as an enhancement, and is so treated in TempLS V4. So the agreement between enhanced but otherwise inferior methods, and the better methods that are not enhanced, is further corroboration of the convergence of integration methodology.

I will illustrate all this with an active graph, of the kind I have been using for the index results themselves. You can, with mouse, change scales, translate axes, choose various subsets to plot, and also smoothe. An added facility here is that you can switch to difference plots in any combination you choose. I will at some stage post a facility like that for WebGL for doing this.

The post will be long, so I start with a table of contents.

Table of contents

Brief summary of methods

Some recent articles are here, here and here, with links to earlier. The methods I'll list are
  • No weighting - the simple average of stations reporting in each month is used. This is a very poor method, but of interest here because it becomes quite serviceable after SH enhancement.
  • Simple grid. This is the traditional method, where the temperature anomalies within each cell of, say, a lat/lon 2° grid, are averaged, and the global is then the area-weighted average of the cells that have results. Area-weighting relates to the shrinking of area of cells near the poles. It is used by HADCRUT; the paper of Cowtan and Way showed that the effect of accounting for cells without data gave an important correction to trends. I do not now use lat/lon grids, but rather a cubed sphere, or other platonic solid based grid. The alternative is usually icosahedron. In each case I use mappings to make the cells of almost uniform area.
  • Grid with infill. This assigns to empty cells a value based on neighbors. Most recently, I do this by solving a Laplace equation, with the known cell values as boundary condition. That sounds complex, but the simple Southwell relaxation method, , which initially guesses the unknown cells, and then replaces with the average of neighbors until convergence, is quite adequate.
  • Irregular triangular mesh. This has been my workhorse; it is basically finite element quadrature, with linear interpolation within triangles with stations as vertices. I have thought of it as my best method.
  • First order LOESS. This sets up a regular array of nodes (icosahedral), and assigns values to them based on local first order weighted regression with a typically 20 of the nearest stations. The regular array is then simply averaged. I think this is a rival for the best method.

Active plot of results

The idea of active plotting is described here, and my regular example is the monthly plot of indices here. The active plot for this post is here, with details below:


Thursday, April 18, 2019

Description and code of TempLS V4 global temperature program.

It's now time to post the code for TempLS V4. I know that it is a day when many folks who really want to get their message read put it out to the world. I didn't want to steal anyone's thunder; the week just didn't have enough days. I've been writing background articles leading up to the release of V4 of TempLS; here are the links:
ideas about spatial integration
icosahedral grid with equal area mapping
documentation system
details of new methods
tests
math methods

TempLS is a code written in R, dating back to about 2010 - there is a summary of its history up to V3 here. The release post (Aug 2010) for V2 is here; Ver 3 was rolled out in three parts, of which the last, here, links to the earlier parts.

Release of V4 is, as I mentioned in the documentation post, complicated by the fact that I use extensively a library of R functions that I keep, and will need to include. But I hope the documentation will also help; I'll be relying on it for this release.

The release is in three zip files:
  • There is a file of code and some internal data - LScode.zip. To run it, you should extract it in an empty directory; it will create six subdirectories (x,g,i,l,m,n). The R file to run it is LS_wt.r, explained below. This zip file is about 98 Kb.
  • As mentioned before, the code now has a lot of generic functions embedded from my general library. There is a documentation system explained at that link, which I also used for the functions of TempLS, below. This zip file includes a list Moyhupack.rda. This you can attach, R style, as a package, which has advantages. The functions from it are also on a file Moyhupack.r, which you can just source() into your environment. It has a lot of functions with short names, most of which you won't want to know about, so there is potential for naming clashes. Finally, there is a documentation file Moyhupack.html. If there is interest, I will write more about the library. I think it does have useful things for math programming. The zipfile is about 800 Kb
  • Finally, there is a set of data. The main one is a recent version of GHCN V4 unadjusted TAVG, named ghcnv.dat. There is also the inventory file that came with it, called ghcnv.inv. The SST data is from ERSST V5; it is part processed in TempLS style and in a directory called x. Note that the code zipfile also created an x directory, with just one list, which would overwrite this larger version if unzipped later. This is a big file - about 54 Mb (mainly GHCN).

The math task and code structure

As before (V3 and earlier), the code is divided into four parts, The code starts with a text version of GHCN V4, and goes looking for a NetCDF version of ERSST V5, if needed. It assembles all this data into 12 (month) arrays of station temperature (rows) vs year (cols). ERSST grid cells count as stations. Missing values are marked with NA. The objective is to fit the model
T=L+G
where T is the station temperature, L the time-invariant normal, and G the space-invariant global anomaly. The fitting requires integration, which comes down to a weighted sum of the data. The third section of the code uses the data about stations reporting, but not the temperatures, to work out weights for integration by various methods. This can be by far the most time-consuming part of the code, alleviated greatly by the use of past calculated weights where the list of stations reporting hasn't changed.

The fourth section is where the fitting is done, via a brief iteration (few seconds) to get convergent values of L and G, which is the main output. It is also where the optional Spherical Harmonics enhancement is done.

Code details.

The code is now almost entirely done as function calls. I call the main file LS_wt.r, which calls functions from LS_fns.r. These are documented below. The main sequence functions (the four parts) are
  • SortLandData() handles land data
  • SortSSTData() handles ERSST5 data.
  • DeriveWeights() calculates weights from the x array of the previous parts
  • FitGlobal() does the iterative fitting to produce L and G.
There are three lists of data that are switched in and out as needed; they are:
  • I the inventory data, which governs the dimensions of later data matrices. It is updated in a controlled way by newinv() on LSinv.rds

  • J is to describe the results in SortLandData and SortSSTData - assembling x. On LS.rds
  • K is for the method-dependent data in DeriveWeights and FitGlobal, including the main result, K$G.
The 12 (month) temperature files are stored in directory x. The 12 weight files are stored in directories, one for each method
("g","i","m","n","l") for ("grid","infilled","mesh","none","loess"). Since the K data is method dependent, a separate version is stored on each directory as LS.rds.

Code sequence

The main code is:
source("LS_fns.r")
if(!exists("job"))job="UEMO"
print(paste("Starting",job,Sys.time()))
i=vsplit(job,"")
kind=tolower(i)
kinds=c("g","i","m","n","l")
RDS=pc(kind[3],"/LS.rds")
K=readRDS(RDS) # I is inv list, J is list about x, K is list for method kind (w)
K$edition=I$edition
K$kind=kind; K$job=job;
do=i!=kind
wix=pc(kind[3],"/") # info about w
tm=timex(0); tx=unlist(tm); t0=Sys.time();
yr=c(1900,tx[1]+(tx[2]>1)-1); yrs=yr[1]:yr[2]; # you can change these
ny=diff(yr)+1
saveR(K)

if(do[1])SortLandData()

if(do[2])SortSSTData()

if(do[3])DeriveWeights()

if(do[4])K=FitGlobal()


As previously, the user supplies a job string of four characters. They are uppercase if that code section is to be performed. A typical value is job="UEMO". At the moment there aren't realistic alternatives for "UE", which is unadjusted GHCN V4 and ERSST V5. M stands for mesh method, but could be any of ("G","I","M","N","L"). "O" just means perform last section; "P" means go on to produce graphics with another program.

A second control variable that can be set is nSH=4, say. It induces the spherical harmonics enhancement, and the value sets the order, using (nSH+1)^2 functions. Going past 10 or so is risky for stability (and time-consuming).

The third control is called recalc. It lets you override the system of using stored values of weights when the stations reporting in a given month is unchanged. This saves time, but it can happen that you suspect the stored data is wrong, for some reason. Something might have gone wrong, or the inventory might have changed. The default setting is FALSE, or 0, but if you want it to not use stored data, set recalc to 1. There is also a setting that I find useful, recalc=2, which recalculates only the most recent year. This takes very little time, but will check if the code is currently working for that method option. Otherwise if it uses entirely stored data, it could take some time to find errors.

So the actual code here just brings in K, which also acts as a record of what was done. It stores some information and puts K back on disk. The other stuff just makes some time information for use in the main sequence. Note that the last step outputs K. This is where the results are (K$G and K$L).

Documentation of functions

Remember there are a lot of generic functions on the Moyhu package. The functions here are those specific to TempLS.

Tuesday, April 16, 2019

GISS March global up 0.21°C from February.

The GISS V3 land/ocean temperature anomaly rose 0.21°C in March. The anomaly average was 1.11°C, up from December 0.90°C. It compared with a 0.208°C rise in TempLS V3 mesh. Jim Hansen's detailed report is here. So far, April is looking warm too.

I think that now that TempLS and GISS are using GHCN V4, the agreement will be even better than in the past, as in this month. There extra coverage does make a difference. The earlier NCEP/NCAR average also agreed very well (0.19° rise). It was the third highest March in the record, just behind 2017.

The overall pattern was similar to that in TempLS. Huge warm spots in Siberia and NW N America/Arctic. Cool spots in NE USA through Labrador to Greenland, and Arabia through N India. Warm in Australia and S Africa.

As usual here, I will compare the GISS and previous TempLS plots below the jump.

Monday, April 15, 2019

The math basis for the TempLS V4 global temperature program.

This is, I hope, the last of the preparatory posts before I post the code and description of TempLS V4. Earlier posts were on ideas about spatial integration, an icosahedral grid with equal area mapping for a new method, the documentation system that I'll be using, details of new methods, tests, and some math methods that I'll be using and referring to.

The math basis is similar to what I described back in 2010 for V2, and in the guide to V3, still the most complete reference. The statistical model for temperature for a given month (Jan-Dec) is:

T = L + G, where T is measured temperature, L is a time-invariant offset for stations, and G a space-invariant global average temperature. In more detail, using my index notation, with s for station, y for year, and m for month, it is

Tsmy = IyLsm + IsGmy

The I's are added for consistency of indexing; they denote arrays of 1's.

A change in V4 is that the analysis is done separately over months, to avoid over large arrays in RAM. That simplifies the analysis; the subscript m can be dropped. Another change is more subtle. The model fitting should be by area and time, ie by minimising the integral of squares

In earlier versions I discretised the variables jointly, giving a weighted sum

w(sy)(Tsy-IyLs + IsGy)2

In fitting, this gave a matrix which was symmetric positive definite, which seemed natural and good. But it means that if you choose the weights to be right for the spatial integration, they can't be controlled for the time integration, so locally, they vary over years, which isn't really wanted. For the methods I used in v3, the weights were all positive, and of reasonably uniform size, so it didn't really matter. But in V4 I use methods which are higher order, in taking account of derivatives. Sometimes, for good reason, the weights can be negative. Now for spatial integration, they add to the area, so that is still OK. But for time integration, there might be significant cancellation, with a sum that is near zero, or even negative. It doesn't often happen, but since you have to divide by the sum of weights to get a time average, it is bad.

Revised weighting - Gauss Seidel iteration

The previous equations, derived from the sum squares, were
w(s)y(Tsy-IyLs - IsGy)=0(1a)
ws(y)(Tsy-IyLs - IsGy)=0(2a)

To get even time weighting in (1a), I use a weighting J(s)y, which is 1 where there is a reading, and 0 where there isn't (where w also had zeroes).
J(s)y(Tsy-IyLs - IsGy)=0(1b)
ws(y)(Tsy-IyLs - IsGy)=0(2b)
This can be written in part-solved form as
J(s)yIyLs = J(s)y(Tsy - IsGy)=0(1)
ws(y)IsGy = ws(y)(Tsy-IyLs)=0(2)

The first just says that L is the time average of T-G where observed, since its multiplier is just the set of row sums of J, which is the number of years of data for each station. The second then is just the area average of T corrected for L (the anomaly, whne L has converged). Suitably normalised, the multiplier of G is 1. Starting with G=0, this gives "naive averaging", with L just the overall time mean at station s. As I wrote here, this gives a bad result, which is why the main methods insist on averaging over a fixed time period (eg 1961-90). But that then leaves the problem of stations that do not have data there, which this method avoids.

So the first equation is solved, and the L used to estimate G, which is then used to correct, L, and so on iterating. The reason that iteration works is that the equations are weakly coupled. Eq (1) almost fixes L, with variations in G having only a small effect. Conversely Eq (2) almost fixes G, not very sensitive to local variations in L. There is an exception - if you add a constant value to G, it will cause the L's to drop by a similar constant amount.

So that is the iteration sequence, which can be characterised as block Gauss-Seidel. Start with a guessed value for G. Solve (1) for L, then solve (2) for an updated value of G. For GHCNV4, as used here, this converged gaining at least one significant figure per iteration, so four or five steps are sufficient. In practice, I now start with a previous estimate of G, which converges even faster. But in any case, the step takes only a few seconds. At each step, G is normalised to have mean zero between 1961 and 1990 (for each month), to deal with the ambiguity about exchanging a constant value with L.

Program structure

As before, TempLS has four blocks of code, now expressed as functions:
  • SortLandData
  • SortSSTData
  • DeriveWeights
  • FitGlobal
The first two are just housekeeping, ending with xsy, the array of monthly temperatures. The third derives the corresponding spatial integration weight matrix wsy by one of five methods described here. The fourth, FitGlobal(), performs the iteration described above. The result are the parameters G and L, of which G is the desired global temperature anomaly, which I publish every month.

For the more accurate methods, DeriveWeights is the most time consuming step; a full meshing can take an hour, and LOESS takes ten minutes or so to do the full record since 1900. But the weights depend only on the stations reporting, not what they report, and for most of those years this doesn't change. So I store weights and reuse them unless there is evidence of change in the list of stations that reported that year. This brings the compute time back to a few seconds.

In V3, I had a separate method based on Spherical Harmonics. As described here, I now treat this as an enhancement of any method of integration. In V3, it was in effect an enhancement of the very crude method of unweighted averaging over space. It actually worked well. In V4 it is implemented, optionally, at the start of FixGlobal(). The weights from part 3 are modified in a quite general way to implement the enhancement, with results described in the post on tests. I think it is of interest that as a separate integration principle, it yields results very similar to the more accurate (higher order) integration methods. But I don't think it will have a place in routine usage. It takes time, although I now have a much faster method, and it does not give much benefit if the more accurate methods are used. So why not just use them?







Friday, April 12, 2019

Some math used in TempLS - index notation and sparse matrices.

Index notation for arrays, and the summation convention

Index notation for arrays became popular with the tensor calculus, partly because it also elegantly embraced the concepts of contravariance and covariance. It is sometimes called the Einstein notation. I'm not using the tensor aspect here, so no superscripts. Arrays, including vectors, are simply represented by writing down their coefficients with subscripts, with the understanding that those vary over a known range. There is no attempt to visualise them geometrically. The coefficients, like arrays, can be added if the subscripts match in range, and be multiplied by scalars or pother subscripted arrays.

But the core of the system is the summation convention, also named after Einstein. It expresses the idea of inner product, which is what distinguishes matrices from just vectors. If an index is repeated in a term, it implies summation over the range of that index. Some familiar cases:

aibjis the outer product of two vectors, functioning as a 2-index array, or matrix, but
aibiis the inner product, with a single number (scalar) result. i is referred to as a dummy variable, because it won't be referenced again.
aijbjis a matrix right multiplied by a vector. The result is a vector, indexed with i
aijbiwould be left multiplied. But that terminology isn't needed; the indices say it all.
ajjis the trace of matrix A
aijbjkis the matrix product A*B, indices i and k.
aijbjkckis the product A*B*c, a vector, index i


You can exempt an index from the repetition count with braces, so a(i)bi is not summed. The multiplication is real, so it is as if A was a diagonal matrix. It often appears as a kernel, as in a(i)bici, which would be a quadratic form with diagonal matrix A as kernel.

I use a special notation Ij or Ijk for an array consisting entirely of 1's.

I have used this notation at least since introducing V2 of TempLS, but I'll be making more extensive use of it with V4. A feature is that once you have written down an index version of what you want to do, it maps directly onto R coding.

Sparse matrices

Often the indices have a small range, as over space dimensions. But sometimes the range is large. Most integration methods in TempLS involve mapping from one large set, like stations, to another, like grid cells. At some point there will be a matrix involving those two indices. In V4 we might have 30,000 stations, and 10,000 cells. You would not want to enter such a matrix into RAM.

The saving grace is that most nodes have no interaction with most cells. As with partial differential equations in maths, the relations are local. Most terms are zero, a sparse matrix. So it is possible to enter the non-zero terms, but you also have to say where they are in the matrix. I just list the N nonzero terms in a vector and in a Nx2 matrix, I list first the row and then the column of each term.

The key issue again is matrix-vector multiplication y=A*b. This is what can take a set of numbers over stations to one over grids. In a language like C, with this structure, the algorithm is simple. You go through that list, take each apq, multiply by bq, and put the result in yp. But that is a lot of operations, and in an interpreted language like R, they carry too much overhead. It would be possible to stretch out the terms of b to match a, so the multiplication could be between vectors. But they can't be added into the result vector in one operation, because several numbers would have to be added into one location. I use an algorithm that sorts the indices of a into subsets that do correspond to a unique location, so the multiplication can be reduced to a small number of vector operations.

Use in TempLS

I have written. for example here and very recently here about ways in which spatial integration is done in TempLS, and why. Integration is a linear operation on a set Ts of measured temperatures, or anomalies, and so the result can be represented as a weighted sum, wsTs. For each method, the first task of TempLS is to calculate those weights.

I'll go through the basic arithmetic of simple gridding from this viewpoint. To integrate directly, you work out which cells of the grid the stations belong to, calculated an average for each cell, and then multiply those by the areas of each cell and add. It's fairly simple to program, but I'll describe it in sparse matrix terms, because more complex cases have a similar structure.

The key matrix being formed is the incidence matrix In over indices g for grid (rows) and s for stations (columns). This is zero, except where station s is in cell g, when it has a 1. So the action of summing station data x in cells is the product Ingsxs. The sum has to be divided by the count, which in matrix terms is the product IngsIs. The result then has to be multiplied by the grid areas ag and added. Overall

Integral over earth = ag (ItIn(g)t)-1 Ings xs

The result has no indices, so all must be paired. Note the appearance of a bracketed (g) in the denominator. The repetition of separate station indices s and t indicates two separate summations.

Forming the index set for the sparse matrix is usually simple. Here the right column just numbers the stations, starting from 1. The left column just contains the cell number that that station is in.

That was reasoned starting from known x. To calculate the weights w, x must be removed, and the linear algebra started from the other end, with matrices transposed. Transposing a sparse matrix is trivial - just swap the columns of the index.

I'll go through the algebra for more complex cases in a future post, probably in conjunction with the code. But here is a table showing how the basic structure repeats. These are the various integration algorithms in index notation.

No weighting (ItIt)-1Isxs
Simple grid ag (ItIn(g)t)-1 Ings xs
Triangle mesh ag Ings xsa = areas of triangles
Loess order 0 Ip (ItW(p)t)-1 Wps xsW is a function (with cutoff) of distance between station s and nearby node p.
Loess order 1 ypj (ztjztkW(p)t)-1 Wpszsk xs zsk are the station of station s; ypj are the coordinates of node p.






Tuesday, April 9, 2019


Testing integration methods in TempLS V4.

I posted yesterday about some new methods and changes in integration in TempLS V4, explaining how the task is central to the formation of global temperature anomaly averages. Today I will show the results of testing the methods. We can't test to see if actual data gives the right result, because we don't know that result. So tests are of two types
  • Testing integration of functions whose average is known. Values of those functions replace the observed temperatures in the process, with the same station locations and missing values.
  • Comparing different methods, to see if some converge on a common answer, and which do not.

I tested six methods. They are, in order (approx) of ascending merit
  • Simple average of stations with no area weighting.
  • Simple gridding, with cells without data simply omitted. The grid is cubed sphere with faces divided into 24x24, with a total of 3456 cells. Cell area is about 25000 sq km, or 275 km edge.
  • Zero order LOESS method, as described in last post. Zero order is just a local weighted average.
  • Gridding with infill. Empty cells acquire data from neighboring stations by solving a diffusion equation, as briefly described in the last post.
  • First order LOESS, in which values on a grid of nodes are calculated by local linear regression
  • Finite element style integration on a triangular mesh with observation points at the vertices.
I said the order was approximate; there is some indication that full LOESS may be even slightly better than mesh. I also tested the effect of using spherical harmonics fits for enhancement, as described here. This option is now available in TempLS V4. The parameter nSH is a measure of the number of periods around the Earth - see here for pictures. For each level of nSH, there are (nSH+1)^2 functions.

Testing known functions

An intuitively appealing is simply latitude. Using just the latitude of stations when they report, in January 2011, what does the process give for the average latitude of the Earth. It should of course be zero.


nSH=0nSH=2nSH=4nSH=8nSH=12
No weighting26.6295-0.0393-0.06030.04791.9308
Grid no infill0.50950.0806-0.0121-0.035-0.0197
LOESS order 00.022-0.0342-0.0438-0.0342-0.0216
Grid diffuse-0.0432-0.0466-0.037-0.0267-0.0175
LOESS order 1-0.023-0.0231-0.0252-0.0223-0.0182
Mesh FEM-0.0209-0.022-0.0218-0.0177-0.0129


The case with no weighting gives a very bad result. A very large number of GHCN V4 stations are in the US, between lat 30° and 49°, and that pulls the average right up. Grid with no infill also errs on the N side. The reason here is that there are many cells in the Antarctic and near without data. Omitting them treats them as if they had average latitude (about 0), whereas of course they should be large negative. The effect is to return a positive average. The other methods work well, because they are not subject to this bias. They have errors of local interpolation, which add to very little. The poor methods improve spectacularly with spherical harmonic (SH) enhancement. This does not fix the biased sampling, but it corrects the latitude disparity that interacts with it. Deducting the non-constant spherical harmonics, which are known to have average zero, leaves a residual with is not particularly biased by hemisphere. The no weighting case is bad again at high order SH. The reason is that the fitting is done with that method if integration, which becomes unreliable for higher order functions. I'll say more about that in the next section.

Testing spherical harmonics.

Latitude is a limited test. SH's offer a set of functions with far more modes of variation, but testing them all would return a confusing scatter of results. There is one summary statistic which I calculate that I think is very helpful. When you fit SH Fourier style, it is actually a regression, involving the inverse of the matrix of integrals of products of SH. This is supposed to be the identity matrix, but because of approximate integrations, the SH are not exactly orthogonal. The important number indicating deviation from unit is the condition number, or ratio of max to min eigenvalues. When it gets too low, the inversion fails, as does the fitting. That is what was happening to unweighted averaging with nSH=12. So I calculated the minimum eigenvalue of that matrix (max remains at 1). Here is a table:

nSH=2nSH=4nSH=8nSH=12
No weighting0.07030.0205-0.8361-0.9992
Grid no infill0.78930.6310.43540.1888
LOESS order 00.95010.82980.45930.1834
Grid diffuse0.96940.90130.60260.2588
LOESS order 10.9880.95030.69050.3037
Mesh FEM0.98750.94180.66770.3204


Close to 1 is good. There is a fairly clear ranking of merit, with full LOESS just pipping mesh. No weighting is never good, and completely falls apart with high nSH, being not even positive definite. All methods are down quite a lot at nSH=12, although a minimum of 0.3 say is not a major problem for inversion. These values are for January 2011, but it is a robust statistic, varying little from year to year. Here is the comparison over years for the nSH=8 level:


201120122013201420152016201720182019
No weighting-0.8361-0.759-0.2471-0.2119-0.0423-0.01040.010.01670.0199
Grid no infill0.43540.4650.44810.41290.44930.42370.43940.43250.4436
LOESS order 00.45930.49720.47420.44770.47910.44720.47340.46970.4764
Grid diffuse0.60260.62070.61020.57410.61320.59290.61270.6080.6095
LOESS order 10.69050.7110.69650.67620.69120.68320.69030.69270.6944
Mesh FEM0.66770.66630.67030.66360.66270.67370.67410.67720.6626


Test of global average of real anomalies.

For consistency, I have used real temperature data with the station normals calculated by the mesh method. It would not matter which method was used, as long as it is consistent. Here are January results over years with no SH enhancement:


201120122013201420152016201720182019
No weighting-0.15651.22770.53410.34751.04010.72130.98810.80650.8183
Grid no infill0.32540.29950.48340.52510.67290.94640.81940.61270.7443
LOESS order 00.3660.29880.52240.56650.66491.00360.84990.6190.7425
Grid diffuse0.36740.29190.51130.54430.64950.98020.82830.61860.7111
LOESS order 10.37060.29940.50570.53970.64290.99320.82810.62090.7112
Mesh FEM0.37790.28620.50510.5170.64740.99690.80880.61780.7017


The best methods are really very close, generally within a range of about 0.02. Even simple grid (with cubed sphere) is not that different. But no weighting is bad. The effects of SH enhancement are variable. I'll show them for January 2018:


nSH=0nSH=2nSH=4nSH=8nSH=12
No weighting0.80650.88420.56870.56331.4847
Grid no infill0.61270.62530.60770.61650.6139
LOESS order 00.6190.61960.61380.62240.6152
Grid diffuse0.61860.61940.61910.62410.6237
LOESS order 10.62090.6210.61980.62250.619
Mesh FEM0.61780.61820.6180.62060.6147


The better methods do not change much, but converge a little more, and are joined by the lesser methods. This overall convergence based on the separate principles of discretisation type (mesh, grid etc) and SH enhancement is very encouraging. Even no weighting becomes respectable up to nSH=8, but then falls away again as the high order fitting fails. I'll show in a future post the comparison results of the different methods for the whole time series. There are some earlier comparisons here, which ran very well back to 1957, but were then troubled by lack of Antarctic data. I think LOESS will perform well here.

Monday, April 8, 2019


New methods of integration in TempLS V4 for global temperature.


Background

TempLS is a program that takes the extensive data of surface temperature measurements and derives a global average of temperature anomaly for each month over time. It also produces maps of temperature anomaly distribution. The basic operation that enables this is spatial integration. As with so much in science and life, for Earth's temperature we rely on samples - it's all we have. To get the whole Earth picture, it is necessary to interpolate between samples. Integration is the process for doing that, and then adding up all the results. The average is the integral divided by the area.

The worst way of getting an average is just to add all the station results and divide by the total. It's bad because the stations are unevenly distributed, so the result reflects the regions where stations are dense. This generally means the USA. Some kind of area weighting is needed so that large areas with sparse readings are properly represented. Early versions of TempLS used the common method of gridding based on latitude/longitude. The default method of spatial integration is to form a function which can be integrated, and which conforms in some way to the data. In gridding, that function is constant within each cell, and equal to the average of the cell data. But there is the problem of cells with no data...

Since V2, 2011 at least, TempLS has used unstructured mesh as its favored procedure. It is basically finite element integration. The mesh is the convex hull of the measurement points in space, and the area weight is just the area of triangles contacting each node. For over seven years now I have reported average temperature based on the mesh method (preferred) and grid, for compatibility with Hadcrut and NOAA.

Early in the life of V3, some new methods were added, discussed here. The problem of cells with missing data can be solved in various ways - I used a somewhat ad-hoc but effective method I called diffusion. It works best with grids that are better than lat/lon. I also used a method based on spherical harmonics, with least squares fitting. As described here, I now think this should be seen as an enhancement which can be applied to any method. It is spectacularly effective with the otherwise poor method of simple averaging; with better methods like mesh or diffusion, there is much less room to improve.

So why look for new methods?

We don't have a quantitative test for how good a method is applied to a temperature field. The best confirmation is that the methods are relatively stable as parameters (eg grid size) are varied, and that they agree with each other. We have two fairly independent methods, or three if you count SH enhancement. They do agree well, but it would be good to have another as well.

V4 changes.

V4 does introduce a new method, which I will describe. But first some more mundane changes:
  • Grid - V4 no longer uses lat/lon grids, but rather grids based on platonic solids. Currently most used is the cubed sphere, with ambitions to use hexagons. All these grids work very well.
  • Spherical Harmonics - is no longer a separate method, but an enhancement available for any method. It's good to have, but adds computer time, and since it doesn't much enhance the better methods, it can be better to use them directly.
  • I have upgraded the diffusion method so that it now solves a diffusion equation (partial differential) for the regions without cell data. The process is very simple - Southwell relaxation, from the pen and paper era, when computer was a job title. You iterate replacing unknown values by an average of neighbors.

The LOESS method

The new method uses local regression - the basis of LOESS smoothing. Other descriptive words might be meshless methods and radial basis functions. The idea is that instead of integrating the irregular pattern of stations, you find a set of regularly spaced points that can be integrated. In fact, using an icosahedron, you can find points so evenly spaced that the integral is just a simple average. To estimate the temperatures at these points, weighted regression is applied to a set of nearby measurements. The regression is weighted by closeness; I use an exponential decay based on Hansen's 1200 km for loss of correlation. But I also restrict to the 20 closest points, usually well within that limit.

The regression can be relative to a constant (weighted mean) or linear. The downside of constant is that there may be a trend, and the sample points might cluster on one side of the trend, giving a biased result. Linear fitting counters that bias.

I'll show test results in the nest post. I think the LOESS method is at least as accurate as the mesh method, which is to say, very accurate indeed. And of course, it agrees well. It is flexible, in that where data is sparse, it just accepts data from further afield, which is the best that can be done. You could think of a grid method as similarly estimating the central values, which can then be integrated. The grid method, though, artificially cuts off data that it wall accept at the cell boundary.

The LOESS method also gives a good alternative method of visualisation. My preferred WebGL requires triangles with values supplied at corners, when GL will shade the interior accordingly. I have used that with convex hull mesh (eg here), but when triangles get large, it produces some artifices. Using the underlying icosahedral mesh of LOESS has uniformly sized triangles. Of course, this is in a way smoothing over the problem of sparse data. But at least it does it in the best possible way.

Here is a WebGL plot of February 2019 temperature anomaly, done the LOESS way. As usual, there are checkboxes you can use to hide the mesh overlay, or the colors, or even the map. More on the facility and its use here.



You can contrast the effect of the LOESS smoothing with the unstructured mesh representation here. Both present unadjusted GHCN V4, which clearly has a lot of noise, especially in the USA, where quantity seems to degrade quality. None of this detracts from global integration, which smooths far more than even LOESS. I think that which it is occasionally of interest to see the detail with the mesh, the LOESS plot is more informative. The detail of mesh had been useful in GHCN V3 for spotting irregularities, but in the US at least, they are so common that the utility fades. In much of the rest of the world, even Canada, coherence is much better.