Hej Erik,

A report on a bug in the GribAPI, discovered at SMHI. Could you please
forward this to the relevant person. Further contacts concerning this
issue can be directly with our developer Magnus Berg (magnus.berg@2xper.se).

Ha det bra
Ilmar

The function grib_nearest_distance in grib_nearest.c calculates
incorrect distances between the two points. The end result being that
the 4 nearest points returned from grib_nearest_find gets incorrect
distances.

The debug printout below shows the 4 points returned when asked for the
point (-90,0.24). As you can see, the distance from that point to both
points (-89.75,0.25) and (-89.75,0.0) are identical eventhough it should
be closer to (-89.75,0.25), but the worst is that the distances to both
(-90,0.25) and (-90,0.0) are 0. When selecting the nearest point from
these 4, you will get either the 3rd or 4th depending on the algorithm
used to select one, also, if you want to interpolate a value for the
actual point from the four nearest points, the interpolated value will
be incorrect.

distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-89.750000, lon2 0.250000, a 0.999990, a2 0.999990, ret 27.783329
distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-89.750000, lon2 0.000000, a 0.999990, a2 0.999990, ret 27.783329
distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-90.000000, lon2 0.250000, a 1.000000, a2 1.000000, ret 0.000000
distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-90.000000, lon2 0.000000, a 1.000000, a2 1.000000, ret 0.000000

The error is the calculation of a in the function
(grib_nearest_distance), it says:

a = sin(rlat1)*sin(rlat2)+ cos(rlat1)*cos(rlat2)*cos(rlon2-rlon1);

but should be:

a = cos(rlat1)*cos(rlat2) + sin(rlat1)*sin(rlat2)*cos(rlon2-rlon1);

Then the distances become correct, see:

distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-89.750000, lon2 0.250000, a 0.999990, a2 0.999990, ret 27.805547
distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-89.750000, lon2 0.000000, a 0.999982, a2 0.999982, ret 38.513689
distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-90.000000, lon2 0.250000, a 1.000000, a2 1.000000, ret 1.111333
distance: radius 6367.470000, lat1 -90.000000, lon1 0.240000, lat2
-90.000000, lon2 0.000000, a 0.999991, a2 0.999991, ret 26.671996

The bug is present in all the versions I've checked, 1.9.8, 1.9.9 and
1.9.16.

After some further testing, I also see a problem around the equator, it
might be a sign problem when the latitude sign changes. Let me know if
you want more info about that or if you can reproduce it on your side.

As a side note, is it really better to use sphere calculations to
calculate these distances? Wouldn't it be just as good to use linear
trigonometry (Pythagoras)? It's much faster and probably not
significantly different. The earth isn't spherical anyway so the radius
is different on different latitudes.

Regards,

Magnus Berg