Haversine formula returning NaN

[Kevin Mess]

For you mathematicians out there :)

I have a class whose primary attributes are two doubles representing a
geographic point, expressed as a latitude and longitude.

One of the methods in the class calculates the distance between two
Points using the Haversine Formula.

For the most part, it works, but in some instances, I get NaN (Not a
Number) where I realistically should expect a distance.

Here is the method that implements the Haversine formula:

public double distanceTo(Point p2) {
/*
Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
Sky and Telescope, vol. 68, no. 2, 1984, p. 159):
http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html

dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin^2(dlat/2) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c

The Haversine Formula can be expressed in terms of a two-argument
inverse tangent function instead of an inverse sine. No bulletproofing is
needed for an inverse tangent (theoretically at least).
*/

final double RADIUS = 6371.005076123; // Earth's mean radius in km

// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);

// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);

double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) *
sin(dLon / 2);

double c = 2 * atan2(sqrt(a), sqrt(1 - a)); // great circle
distance in radians
return RADIUS * c;
// great circle distance in km
}

I won't provide you with the entire output for brevity. Notice the
first few lines - they return expected results. The last few show the
problem:

Distance from Las Vegas (N36°10'59" W115°13'00") to:
// these work
offshore Northern California (N40°18'29" W124°22'49") = 401.5591308690782 km.
Nias region, Indonesia (N00°47'03" E097°19'49") = 4667.779825988913 km.
offshore Northern California (N40°18'35" W124°23'04") = 401.64873756773324 km.
Gulf of Alaska (N58°20'18" W148°17'06") = 2430.6375616198306 km.
// these don't work
Oklahoma (N34°27'52" W097°47'44") = NaN km.
eastern Turkey (N39°30'18" E040°42'10") = NaN km.
Central California (N35°51'19" W120°24'25") = NaN km.
Caspian Sea, offshore Turkmenistan (N40°25'44" E052°01'29") = NaN km.
...

Any thoughts on where/why this formula is failing?

Thanks,
Kevin

[Patricia Shanahan]

Post by Kevin Mess
For you mathematicians out there :)
I have a class whose primary attributes are two doubles representing a
geographic point, expressed as a latitude and longitude.
One of the methods in the class calculates the distance between two
Points using the Haversine Formula.
For the most part, it works, but in some instances, I get NaN (Not a
Number) where I realistically should expect a distance.
public double distanceTo(Point p2) {
/*
Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin^2(dlat/2) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c
The Haversine Formula can be expressed in terms of a two-argument
inverse tangent function instead of an inverse sine. No bulletproofing is
needed for an inverse tangent (theoretically at least).
*/
final double RADIUS = 6371.005076123; // Earth's mean radius in km
// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);
// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);
double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) *
sin(dLon / 2);
double c = 2 * atan2(sqrt(a), sqrt(1 - a)); // great
circle distance in radians
return RADIUS *
c; // great circle
distance in km
}

...

Post by Kevin Mess
Any thoughts on where/why this formula is failing?

If a is negative, then Math.sqrt(a) is a NaN. Similarly, sqrt(1-a) is a
NaN if a is greater than 1.

If either input to atan2 is a NaN, so is the result, c, and RADIUS*c.

What constrains a to the range 0 through 1?

In the formula in the comment, sqrt(1-a) was never calculated, and there
was a min(1,sqrt(a)) which makes no sense if a cannot be greater than 1.
The square root of a number in the range 0 through 1 is no greater than 1.

Patricia

[Fred Kleinschmidt]

Post by Patricia Shanahan

Post by Kevin Mess
For you mathematicians out there :)
I have a class whose primary attributes are two doubles representing a
geographic point, expressed as a latitude and longitude.
One of the methods in the class calculates the distance between two
Points using the Haversine Formula.
For the most part, it works, but in some instances, I get NaN (Not a
Number) where I realistically should expect a distance.
public double distanceTo(Point p2) {
/*
Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin^2(dlat/2) + cos(lat1) * cos(lat2) *
sin^2(dlon/2)
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c
The Haversine Formula can be expressed in terms of a two-argument
inverse tangent function instead of an inverse sine. No bulletproofing is
needed for an inverse tangent (theoretically at least).
*/
final double RADIUS = 6371.005076123; // Earth's mean radius in km
// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);
// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);
double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) *
sin(dLon / 2);
double c = 2 * atan2(sqrt(a), sqrt(1 - a)); // great
circle distance in radians
return RADIUS * c;
// great circle distance in km
}

...

Post by Kevin Mess
Any thoughts on where/why this formula is failing?

If a is negative, then Math.sqrt(a) is a NaN. Similarly, sqrt(1-a) is a
NaN if a is greater than 1.

If computed correctly, 'a' cannot be negative.
Note that

a = sin(x)*sin(x) + cos(y)*cos(z) *sin(w)*sin(w)

The First cannot be negative (it's a squared number).
Neither cosine can be negative, since latitudes y and z
are in [-90, 90] degrees.
So the second term cannot be negative either.

However, This formula does allow 'a' to exceed 1.0

Post by Patricia Shanahan
If either input to atan2 is a NaN, so is the result, c, and RADIUS*c.
What constrains a to the range 0 through 1?
In the formula in the comment, sqrt(1-a) was never calculated, and there
was a min(1,sqrt(a)) which makes no sense if a cannot be greater than 1.
The square root of a number in the range 0 through 1 is no greater than 1.
Patricia

[Patricia Shanahan]

Post by Fred Kleinschmidt

Post by Patricia Shanahan

Post by Kevin Mess
For you mathematicians out there :)
I have a class whose primary attributes are two doubles representing a
geographic point, expressed as a latitude and longitude.
One of the methods in the class calculates the distance between two
Points using the Haversine Formula.
For the most part, it works, but in some instances, I get NaN (Not a
Number) where I realistically should expect a distance.
public double distanceTo(Point p2) {
/*
Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin^2(dlat/2) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c
The Haversine Formula can be expressed in terms of a two-argument
inverse tangent function instead of an inverse sine. No bulletproofing is
needed for an inverse tangent (theoretically at least).
*/
final double RADIUS = 6371.005076123; // Earth's mean radius in km
// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);
// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);
double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) *
sin(dLon / 2);
double c = 2 * atan2(sqrt(a), sqrt(1 - a)); // great
circle distance in radians
return RADIUS * c;
// great circle distance in km
}

...

Post by Kevin Mess
Any thoughts on where/why this formula is failing?

If a is negative, then Math.sqrt(a) is a NaN. Similarly, sqrt(1-a) is a
NaN if a is greater than 1.

If computed correctly, 'a' cannot be negative.

That "if computed correctly" is important because of Eric Sosman's
observation that the cosine arguments this.latitude and p2.latitude are
not converted to radians. Presumably, they are in degrees with numeric
range -90 through +90, and there are certainly angles in the range -90
radians through +90 radians that have negative cosines.

Post by Fred Kleinschmidt
However, This formula does allow 'a' to exceed 1.0

I'm not sure of that. To get sin^2(dlat/2) and sin^2(dlon/2) close to 1
requires the differences in latitude and longitude to be close to 180
degrees, so we need something like lat1 close to plus 90 degrees, lat2
close to minus 90 degrees, and similarly for longitude. But that makes
the cosines of the latitudes close to 0.

I have not worked it out, but I suspect that if it were computed
correctly with perfect precision, a greater than 1 may not be possible.

Patricia

[Kevin Mess]

Post by Patricia Shanahan

Post by Fred Kleinschmidt

Post by Patricia Shanahan

Post by Kevin Mess
For you mathematicians out there :)
I have a class whose primary attributes are two doubles representing a
geographic point, expressed as a latitude and longitude.
One of the methods in the class calculates the distance between two
Points using the Haversine Formula.
For the most part, it works, but in some instances, I get NaN (Not a
Number) where I realistically should expect a distance.
public double distanceTo(Point p2) {
/*
Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin^2(dlat/2) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c
The Haversine Formula can be expressed in terms of a two-argument
inverse tangent function instead of an inverse sine. No bulletproofing is
needed for an inverse tangent (theoretically at least).
*/
final double RADIUS = 6371.005076123; // Earth's mean radius in km
// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);
// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);
double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) *
sin(dLon / 2);
double c = 2 * atan2(sqrt(a), sqrt(1 - a)); // great
circle distance in radians
return RADIUS * c; // great circle distance in km
}

...

Post by Kevin Mess
Any thoughts on where/why this formula is failing?

If a is negative, then Math.sqrt(a) is a NaN. Similarly, sqrt(1-a) is a
NaN if a is greater than 1.

If computed correctly, 'a' cannot be negative.

That "if computed correctly" is important because of Eric Sosman's
observation that the cosine arguments this.latitude and p2.latitude are
not converted to radians. Presumably, they are in degrees with numeric
range -90 through +90, and there are certainly angles in the range -90
radians through +90 radians that have negative cosines.

Post by Fred Kleinschmidt
However, This formula does allow 'a' to exceed 1.0

I'm not sure of that. To get sin^2(dlat/2) and sin^2(dlon/2) close to 1
requires the differences in latitude and longitude to be close to 180
degrees, so we need something like lat1 close to plus 90 degrees, lat2
close to minus 90 degrees, and similarly for longitude. But that makes
the cosines of the latitudes close to 0.
I have not worked it out, but I suspect that if it were computed
correctly with perfect precision, a greater than 1 may not be possible.
Patricia

Patricia - I can't tell you how many times I've looked at that code and
didn't realize I failed to convert to radians! That solved the
problem. There's nothing like a fresh set of eyes to see the obvious :)

To everyone else that answered, thank you very much. It's nice to have
so many helpful folks around.

Kevin

[Fred Kleinschmidt]

Post by Patricia Shanahan

Post by Fred Kleinschmidt

Post by Patricia Shanahan

Post by Kevin Mess
For you mathematicians out there :)
I have a class whose primary attributes are two doubles representing a
geographic point, expressed as a latitude and longitude.
One of the methods in the class calculates the distance between two
Points using the Haversine Formula.
For the most part, it works, but in some instances, I get NaN (Not a
Number) where I realistically should expect a distance.
public double distanceTo(Point p2) {
/*
Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin^2(dlat/2) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c
The Haversine Formula can be expressed in terms of a two-argument
inverse tangent function instead of an inverse sine. No bulletproofing is
needed for an inverse tangent (theoretically at least).
*/
final double RADIUS = 6371.005076123; // Earth's mean radius in km
// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);
// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);
double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) *
sin(dLon / 2);
double c = 2 * atan2(sqrt(a), sqrt(1 - a)); // great
circle distance in radians
return RADIUS * c; // great circle distance in km
}

...

Post by Kevin Mess
Any thoughts on where/why this formula is failing?

If a is negative, then Math.sqrt(a) is a NaN. Similarly, sqrt(1-a) is a
NaN if a is greater than 1.

If computed correctly, 'a' cannot be negative.

That "if computed correctly" is important because of Eric Sosman's
observation that the cosine arguments this.latitude and p2.latitude are
not converted to radians. Presumably, they are in degrees with numeric
range -90 through +90, and there are certainly angles in the range -90
radians through +90 radians that have negative cosines.

Post by Fred Kleinschmidt
However, This formula does allow 'a' to exceed 1.0

I'm not sure of that. To get sin^2(dlat/2) and sin^2(dlon/2) close to 1
requires the differences in latitude and longitude to be close to 180
degrees, so we need something like lat1 close to plus 90 degrees, lat2
close to minus 90 degrees, and similarly for longitude. But that makes
the cosines of the latitudes close to 0.
I have not worked it out, but I suspect that if it were computed
correctly with perfect precision, a greater than 1 may not be possible.
Patricia

If sin() and cos() were infinitely precise, all would be OK.
But roundoff/truncation errors make it possible for 'a' to
slightly exceed 1.0. Hence the min(1., sqrt(a))

[Thomas Weidenfeller]

Post by Kevin Mess
Any thoughts on where

Use a debugger to find out. My guess is on the calculation of a.

Post by Kevin Mess
/why this formula is failing?

My guess is that a either happens to become negative or > 1 due to
rounding errors. Which then lets one of the sqrt() methods return NaN,
which then lets atan2() return NaN. I would further guess there is a
reason for the min(1, sqrt(a)) term in the original formula because of
known possible rounding errors effects.

/Thomas

[Eric Sosman]

Post by Kevin Mess
For you mathematicians out there :)

Not guilty, despite what the diploma says. But ...

Post by Kevin Mess
[...]
// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);
// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);

Observe the conversion to radians. I infer that the
original latitude and longitude are expressed not in radians,
but in whatchamacallits (a technical term we non-mathematicians
like to use a lot).

Post by Kevin Mess
double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) *
sin(dLon / 2);

Note that the arguments to cos() have not been converted
from whatchamacallits to radians. Result: The cosines are not
of the angles you intend, but of some other angles altogether.

[Kevin Mess]

Post by Eric Sosman

Post by Kevin Mess
For you mathematicians out there :)

Not guilty, despite what the diploma says. But ...

Post by Kevin Mess
[...]
// delta (difference) in latitudes
double dLat = radians(p2.latitude) - radians(this.latitude);
// delta (difference) in longitudes
double dLon = radians(p2.longitude) - radians(this.longitude);

Observe the conversion to radians. I infer that the
original latitude and longitude are expressed not in radians,
but in whatchamacallits (a technical term we non-mathematicians
like to use a lot).

Post by Kevin Mess
double a = sin(dLat / 2) * sin(dLat / 2) +
cos(this.latitude) * cos(p2.latitude) * sin(dLon / 2) * sin(dLon / 2);

Note that the arguments to cos() have not been converted
from whatchamacallits to radians. Result: The cosines are not
of the angles you intend, but of some other angles altogether.

Eric, that was it exactly. I looked at this code endlessly, and never
caught it. Fresh eyes are quite helpful :)

Kevin