NASA
World Wind java code is available on the internet. While looking at the their SDK2
EGM96.java file, I discovered that the simple two dimensional linear (
bilenear) interpolation code has an error. The original code:
...
double ul = this.gePostOffset(topRow, leftCol);
double ll = this.gePostOffset(bottomRow, leftCol);
double lr = this.gePostOffset(bottomRow, rightCol);
double ur = this.gePostOffset(topRow, rightCol);
double u = (lon - lonLeft) / INTERVAL.degrees;
double v = (latTop - lat) / INTERVAL.degrees;
double pll = (1.0 - u) * (1.0 - v);
double plr = u * (1.0 - v);
double pur = u * v;
double pul = (1.0 - u) * v;
double offset = pll * ll + plr * lr + pur * ur + pul * ul;
...
I prepared a schematic to visualize the algorithm:
As you can see at the end of above code snippet, the offset formula starts with pll*ll. pll is (1.0 - u) * (1.0 - v). The correct multiplier of the ll term has to be v*(1-u).
A quick
sanity check: Consider the case when lat=latBottom and lon=lonLeft, i.e. the lower left corner. In that case we expect the offset to be equal to ll. Plugging in the values we get v=1, 1-v = 0 and u=0, 1-u = 1. If we use the original (wrong) formula we get offset = 1*0*ll + 0*0*lr + 0*1*ur + 1*0*lu = 0 which is clearly wrong. If we use the corrected formula we get offset = 1*1*ll + 0*0*lr + 0*1*ur + 1*0*lu = ll which is the expected result.
The easiest way to fix the code is to swap the topRow and bottomRow indices when calling getPostOffset() function.
I guess the error was not detected since results do not differ too much because, although the points are used in wrong order, they are close to each other nonetheless, so the result might not have looked suspicious.
I tried to log a bug report on their
Jira site, but I could not since I don't have an account. I was able to file a
bug report on GitHub.
Lessons learnt:
- Thoroughly test 3rd party code, even if it is from NASA or Mathworks. Every code is guilty unless proven innocent. Ignore this advice and you will find yourself chasing strange errors for a very long time.
- Separate interpolation code into a function, write lots of unit test for that function to verify that it works correctly.