averages/integrals on the sphere
How does Ferret compute averages and intervals on the surface of the earth?
Computing averages and integrals on the surface of a sphere requires the application of a correction for the size of grid cells as one moves from the equator to the poles. Prior to Ferret v5.90 Ferret employed a cosine(latitude) correction but this computation has been refined in more recent versions.
We wish to compute the surface area of a grid cell on a sphere, centered at longitude x and latitude y, with angular widths dx & dy (measuring all angles in radians). The surface area of the band between two latitudes y+dy/2, y-dy/2 (see http://mathworld.wolfram.com/Zone.html) is
h = R*sin(y+dy/2) - R*sin(y-dy/2)
is the distance in meters along the earth's axis between the north/south edges of the latitude cell. Interestingly, this area depends only on h and not explicitly on y. The fraction of the band between longitudes x1 and x2 is
So the grid cell area is
A = 2*pi*R*h * dx/(2*pi)
Since sin(a+b)-sin(a-b) = 2cos(a)sin(b), we have
h = 2*R * cos(y) * sin(dy/2)
Thus the grid cell area is
A = R^2 * cos(y) * dx * 2*sin(dy/2)
Prior to v5.90, Ferret computed cell areas on the sphere using a "small dy" approximation:
limit as dy->0 of 2*sin(dy/2) -> dy
A_approx = R^2 * cos(y) * dx * dy
The ratio of this approximate area over the true area is
F = A_approx / A = dy / 2*sin(dy/2) = 1 / sinc(dy/2)
with dy <= pi. This ratio always exceeds 1, which implies that areas were always slightly overestimated prior to v5.90. The fractional error, F-1, exceeds 1e-7 (machine precision) at dy = 0.1deg, 1e-4 at dy = 3 deg, 11% at dy=pi/2, and 57% at dy=pi (largest possible dy). As a result, all integrals prior to v5.90 were slightly too extreme. In addition, averages over grids with non-uniform dy tended to slightly overweight coarsely-gridded latitudes compared to finely-gridded ones. Averages over regions with uniform dy were unaffected by the approximation, so these did not change in v5.90.
Contributed by Andrew Wittenberg