Whoops. There's a third bug in that code.

So, I'm sitting on the bus this morning executing CRU's IDL code in my head when I suddenly realized that there's another more subtle bug in the exact same code I was looking at the other day.

Here's the critical loop once more:

 for i=0.0,nel do begin
  x=cos(a)*(xkm/110.0)*(1.0/cos(!pi*pts1(i,0)/180.0))+pts1(i,1)
  x=x<179.9 & x=x>(-179.9)
  y=sin(a)*(xkm/110.0)+pts1(i,0)
  y=y>(-89.9) & y=y<89.9
  catch,error_value 
  ; avoids a bug in IDL that throws out an occasional 
  ;  plot error in virtual window
  if error_value ne 0 then begin
   error_value=0
   i=i+1
   goto,skip_poly
  endif
   
  polyfill,x,y,color=160
  skip_poly:
 endfor

So, it's plotting those little 32-sided polygons on a flat map of the world and it's making the adjustment to the size so that when the polygon is near the top (or bottom) of the world it gets larger to correctly cover the required area.

But what happens if it plots a polygon near the 'edge of the world'. For example, what happens if it plots a polygon at 85 degrees of latitude and 170 degrees of longitude?

First, here's a picture of a polygon plotted at 85 degrees of latitude but well away from the 'edge of the world'.


Now look at the same polygon at 170 degrees of longitude. See the problem? It doesn't wrap around to the other side. Oops. Since the world is a sphere you'd expect the polygon to reappear on the left hand side of this picture showing the area of influence of the meteorological station being plotted.


So some information is lost for data being plotted near the 180 degrees line. Admittedly, that's in the middle of the Pacific Ocean (although it does cut through some land mass). But if there are any ocean temperature measurements at the 'edge of the world' then bits of their data isn't being taken into account.

I wonder what, if any, impact these three bugs have on the output of this program.

PS. There's actually a fourth problem with this code. The number 110.0. It's being used to convert from kilometres to degrees of longitude and latitude. The same number is used for both even though the Earth isn't a perfect sphere.

The code is using a value of 39,600 km for the circumference of the Earth, whereas the mean value is actually 40,041 km. But, hey, what's an error of 1% between friends?