Thursday, December 03, 2009

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?

7 comments:

Francis Turner said...

Could that bounds issue be the cause of the occasional crashes the catch code is supposed to work around?

Francis Turner said...

you may want to ignore that previous comment, I assume the lines about

x=x<179.9 & x=x>(-179.9)
y=y>(-89.9) & y=y<89.9

are doing the bounds checking (though I don't quite understand how this works)

John Graham-Cumming said...

I wondered about that, but couldn't make my test program crash. I tried plotting all over the place, but with IDL 7.1 it worked fine (except for the wrapping issue mentioned).

John Graham-Cumming said...

The syntax looks odd, but the < is the max operator and the > the min operator. The syntax uses them applied to the x and y vectors to cap bad values.

The & is just a statement separator.

John Graham-Cumming said...

So, I managed to reproduce a crash from the polyfill code. If you remove the check on the range of x.

It's possible to get polyfill to go off into la la land and never return, it's also possible to get it to throw the error: Invalid vertex list.

The nice thing is that when it does work it correctly wraps around solving the edge of world problem.

Jormundgard said...
This comment has been removed by the author.
Jormundgard said...
This comment has been removed by the author.