0

I use the following algorithm to calculate the distance between two points but it's producing unreasonable results. Where am I wrong?

private static double distFrom(double latA, double lngA, double latB, double lngB) {
    double pk = 180/3.14169;

    double a1 = latA / pk;
    double a2 = lngA / pk;
    double b1 = latB / pk;
    double b2 = lngB / pk;

    double t1 = Math.cos(a1)*Math.cos(a2)*Math.cos(b1)*Math.cos(b2);
    double t2 = Math.cos(a1)*Math.sin(a2)*Math.cos(b1)*Math.sin(b2);
    double t3 = Math.sin(a1)*Math.sin(b1);
    double tt = Math.acos(t1 + t2 + t3);

    return 3959*tt;
}
EpicPandaForce
  • 79,669
  • 27
  • 256
  • 428
Sandah Aung
  • 6,156
  • 15
  • 56
  • 98
  • I'm assuming that this is distance on a sphere, not in the plane? Could you link to the formula you are implementing? – gilleain Apr 21 '15 at 14:29
  • 1
    Could you please explain in more detail what your output is, what you expect the output to be and what you have already tried to solve it? – Michael Lang Apr 21 '15 at 14:31
  • 2
    Ok, it's the Haversine formula, and you could just use an implementation like this one : http://rosettacode.org/wiki/Haversine_formula#Java – gilleain Apr 21 '15 at 14:38
  • I am implementing a version of great circle formula suitable for computer programming. I have been trying all kinds of formulas. If you have a version of your own, you are invited to suggest. – Sandah Aung Apr 21 '15 at 14:42
  • What on earth is `3959`? – EpicPandaForce Apr 21 '15 at 14:50
  • 3959 is the radius of the earth in miles – Sandah Aung Apr 21 '15 at 14:51
  • 1
    It would be great if you tell us which values you "give" to this method, which results you get and which results you expect. – Tom Apr 21 '15 at 14:53
  • Miles? I can feel the next Mars Climate Orbiter catastrophe coming from a distance... ( http://en.wikipedia.org/wiki/Mars_Climate_Orbiter#Cause_of_failure ) – EpicPandaForce Apr 21 '15 at 14:55
  • 1
    @EpicPandaForce I actually prefer the use of kilometres but in my country, the Imperial system is still used. – Sandah Aung Apr 21 '15 at 15:07
  • If you need high precision then you could use something like [GeographicLib](http://geographiclib.sourceforge.net/) where the distance in meters can be found via: `Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2, GeodesicMask.DISTANCE).s12` – SamYonnou Apr 21 '15 at 15:46

3 Answers3

7

Replacing 3.14169 with Math.PI will improve a little since the first digits of PI are 3.14159

Also, there is alot of algorithm out there. try http://www.geodatasource.com/developers/java that will give you an answer in miles, kilometers and nautic miles.

JFPicard
  • 5,029
  • 3
  • 19
  • 43
  • 2
    if at least one number is double, the whole result is double, so no, it is not true about integer division. – libik Apr 21 '15 at 14:39
  • 1
    same goes for the return.. make it: return 3959.0*tt; – W vd L Apr 21 '15 at 14:39
  • @libik Your right. See http://stackoverflow.com/questions/3144610/java-integer-division-how-do-you-produce-a-double Anyway, can you explain the 3959 ? – JFPicard Apr 21 '15 at 14:43
  • 2
    @WvdL - not needed, variable tt is double, therefore computation and result is double – libik Apr 21 '15 at 14:45
  • @libik - I assumed that what JFPicard initially said would be true.. i stand corrected. :-) – W vd L Apr 21 '15 at 14:49
  • Sorry for the confusion. That's why I've edited my post, removing misleading information. – JFPicard Apr 21 '15 at 14:50
  • 1
    Nonsense! the coeff in the last line has 4 digits, so, using more than 5 PI digits gives nothing! – Gangnus Apr 21 '15 at 14:51
  • 1
    Of course, there IS an error in Pi, but it is less than rounding error in the last line. – Gangnus Apr 21 '15 at 15:28
1

Better use standard letters - L for longitude, B for latitude. The formula for http://en.wikipedia.org/wiki/Great-circle_distance method you have used is strange. According to wiki it will be:

    double B1 = latA / pk;
    double B2 = latB / pk;
    double dL = (lngA-lngB) / pk;

    double t1 = Math.cos(B1)*Math.cos(B2)*Math.cos(dL);
    double t2 = Math.sin(B1)*Math.sin(B2);

    double tt = Math.acos(t1 + t2);

Less computation -> less rounding errors.

But notice, this acos based formula is correct only mathematically, but not computationally - it has low float precision at short distances! On the mentioned page there are better formulae for computation, based on asin for short distances and a universal one, based on atan. Choose one of them.

Gangnus
  • 24,044
  • 16
  • 90
  • 149
  • There is a need to stick to Java coding rules in my company and the calculated distance doesn't have to be 100% accurate. I was having problems with my algorithm because the results are about 10 times larger than expected. We don't worry about mistakes smaller than one mile. The distances are bound to be incorrect because we are also covering hill areas where great circle algorithms cannot deal with. – Sandah Aung Apr 21 '15 at 15:12
  • 1
    @SandahAung If you use algorithms with bad floating precision, counting with 4 meaning digits you can easily get ZERO meaning digits in result. i.e., it will be totally incorrect. If you cannot count precision alongside with the result, and you obviously can't do it, use only good/well defined/ with high floating precision formulas. – Gangnus Apr 21 '15 at 15:26
  • @SandahAung can you give some test cases in which results are 10 times larger than expected? Just trying out a few points your method seems to give fairly reasonable results (maybe 1-2% error at most). – SamYonnou Apr 21 '15 at 16:08
  • @SamYonnou If the constants have 4 digits, 1% mistake is already 100 times larger than it is acceptable. So, the code has to be changed - Sandah Aung is absolutely right in it. – Gangnus Apr 24 '15 at 07:41
  • @SandahAung 1. You can never have precise results working with floats. But you have to set the acceptable mistake - in absolute or relative form. 2. Java coding rules have no contradictions with any set of accepted terms. For example, nobody uses ordinateOfSatellite, but ySatellite instead. Another problem is if your team doesn't know the standard geodetic terms. THAN yes, use what you do understand. But excuse me, I will use the standard way. 3. Mistakes of one mile for what distance? 4. What is the connection between hill areas and great circle? You DONT work with geoid. – Gangnus Apr 24 '15 at 07:48
  • @Gangnus the great-circle distance has up to a 0.5% error, so OP's error is only 4-5 times greater, given the incorrect PI value used that isnt too bad (correcting the PI value yields a maximum error of around 0.57% which is the same as `gov.nasa.worldwind.geom.LatLon.greatCircleDistance(...)` that I compared it to). If OP can't afford the natural error of a great-circle distance then he should use a different distance method. (Note: I used `net.sf.geographiclib.Geodesic.WSG84.Inverse(...)` as the "true" distance for comparison, and %errors are based on a random sample of 1mil pairs of points) – SamYonnou Apr 24 '15 at 13:09
  • @SamYonnou What do you mean by **mistake** of great circle dist? What do you take as ideal you are comparing against? If we measure the ideal distance on the sphere, GCD is absolutely precise in math sense. But the OP code is not. SandahAung mentions some measures IN HILL AREA. Is he measuring the real distance along these hill surface? Or maybe, along the road that goes through the mentioned hills?.... But even such errors in understanding are not so important as badly defined floating algorithm, that can easily create errors in millions of percents. – Gangnus Apr 24 '15 at 13:24
  • @Gangnus I'm using the WGS84 ellipsoid geodesic distance given by GeographicLib as the ideal (close to what is used for GPS) which does not take into account elevation. – SamYonnou Apr 24 '15 at 14:16
  • @SamYonnou But you do understand, that this distance can be far from the length of the road or a line upon a real mountaneous surface ;-), don't you? So, until we know what OP means by real/ideal distance, we can't talk about MODEL mistakes seriously. So, the only mistake we can discuss, is the floating point algorithm mistake. – Gangnus Apr 24 '15 at 15:05
  • @Gangnus that's definitely true. If OP is comparing their results to actual distances that take into account elevation or perhaps even things like road networks then it makes sense how they are getting errors in the 1000% range. My point was that no matter how much you improve the great-circle distance function you're not going to see an improvement of more than a couple of percent (assuming that no kind of overflow/underflow is happening which I did not notice in my test). – SamYonnou Apr 24 '15 at 15:11
  • @Gangnus Looks like an interesting discussion is going on but let me respond to points that you raised. 1. We currently don't have time to define the acceptable mistake. I think it will evolve over time. 2. We recently introduce coding rules in the company. We had to do this because developers were adhering to all sorts of strange coding norms. Our coding practices include rules such as "all method names and variable names should be camel-cased" and "code should be understandable without reading comments". I was involved in the creation of these rules so for the time being I am my own rules. – Sandah Aung Apr 24 '15 at 18:38
  • 3. We are currently developing apps that would help drivers find our fuel stations. Our job is to build mobile apps that would help them find nearest stations and give them an estimate of how far the stations are. We have / plan to have at least 5 stations in each state and no state in my country is more than 100 miles across so looks like we are talking about a mistake of a mile over a distance of no more than 30 miles. – Sandah Aung Apr 24 '15 at 18:48
  • 4. Let's assume that the nearest station to a user is 1 mile away but there is a mountain between the user and the station so he/she would have to travel a distance of 2 miles to reach there but surprisingly drivers in hill areas seem to be ok with that. Drivers in plain areas however seem to be annoyed by errors smaller than 20 meters because there are fuel stations that could provide fuels matching our quality every few yards so we want accuracy there but even here `acos` doesn't seem to a problem. The problem was the way we were testing the program. Not the algorithm. – Sandah Aung Apr 24 '15 at 18:54
0

Here's an implementation of the atan version from http://en.wikipedia.org/wiki/Great-circle_distance ("special case of the Vincenty formula") which is "accurate for all distances":

/**
 * Mean earth radius in Kilometers (KM) as defined in WGS84
 */
public static final double EARTH_RADIUS_KM = 6371.0087714;

/**
 * Mean earth radius in US Miles, converted from Kilometers
 */
public static final double EARTH_RADIUS_MI = EARTH_RADIUS_KM * 0.621371192;

/**
 * Mean earth radius in Nautical Miles, converted from Kilometers
 */
public static final double EARTH_RADIUS_NM = EARTH_RADIUS_KM * 0.539956803;

/**
 * Calculates the distance between two points on a sphere using a "special
 * case of the Vincenty formula" for accuracy. This uses more trigonometric
 * calls than other formulas, but is more accurate for all distances.
 * <p/>
 * 
 * For more details, see <a href=
 * "https://en.wikipedia.org/wiki/Great-circle_distance">https://en.wikipedia.org/wiki/Great-circle_distance</a>.
 * 
 * @param startLat
 *          The starting point's latitude
 * @param startLon
 *          The starting point's longitude
 * @param endLat
 *          The ending point's latitude
 * @param endLon
 *          The ending point's longitude
 * @param radius
 *          The radius of the spherical earth to be used
 * 
 * @return The distance between the points for the given radius
 */
public static double greatCircleDistance( double startLat, double startLon, double endLat, double endLon, double radius )
{
    // convert to radians
    double lonDelta = Math.toRadians( startLon - endLon );
    double lat1 = Math.toRadians( startLat );
    double lat2 = Math.toRadians( endLat );

    // compute repeatedly used values
    double cos1 = Math.cos( lat1 );
    double sin1 = Math.sin( lat1 );
    double cos2 = Math.cos( lat2 );
    double sin2 = Math.sin( lat2 );
    double cosDelta = Math.cos( lonDelta );

    double top1 = cos2 * Math.sin( lonDelta );
    top1 *= top1; // squared

    double top2 = (cos1 * sin2) - (sin1 * cos2 * cosDelta);
    top2 *= top2; // squared

    double top = Math.sqrt( top1 + top2 );
    double bottom = (sin1 * sin2) + (cos1 * cos2 * cosDelta);

    double rad = Math.atan( top / bottom );

    return radius * rad;
}

Choosing different values for the radius of the earth sets the units desired. This will work for any sphere of any radius, so distances on Mars could be calculated using a radius of 3389.5 km as long as you know the lat/lon of the starting and ending points.

Jason Greanya
  • 189
  • 1
  • 4