How to find distance from the latitude and longitu

Question

I have a set of latitudes and longitudes of locations.

• How to find distance from one location in the set to another?
• Is there a formula ?

Answer 1

The Haversine formula assumes a spherical earth. However, the shape of the earh is more complex. An oblate spheroid model will give better results.

If such accuracy is needed, you should better use Vincenty inverse formula. See http://en.wikipedia.org/wiki/Vincenty's_formulae for details. Using it, you can get a 0.5mm accuracy for the spheroid model.

There is no perfect formula, since the real shape of the earth is too complex to be expressed by a formula. Moreover, the shape of earth changes due to climate events (see http://www.nasa.gov/centers/goddard/earthandsun/earthshape.html), and also changes over time due to the rotation of the earth.

You should also note that the method above does not take altitudes into account, and assumes a sea-level oblate spheroid.

Edit 10-Jul-2010: I found out that there are rare situations for which Vincenty inverse formula does not converge to the declared accuracy. A better idea is to use GeographicLib (see http://sourceforge.net/projects/geographiclib/) which is also more accurate.

Answer 2

Apply the Haversine formula to find the distance. See the C# code below to find the distance between 2 coordinates. Better still if you want to say find a list of stores within a certain radius, you could apply a WHERE clause in SQL or a LINQ filter in C# to it.

The formula here is in kilometres, you will have to change the relevant numbers and it will work for miles.

E.g: Convert 6371.392896 to miles.

    DECLARE @radiusInKm AS FLOAT
    DECLARE @lat2Compare AS FLOAT
    DECLARE @long2Compare AS FLOAT
    SET @radiusInKm = 5.000
    SET @lat2Compare = insert_your_lat_to_compare_here
    SET @long2Compare = insert_you_long_to_compare_here

    SELECT * FROM insert_your_table_here WITH(NOLOCK)
    WHERE (6371.392896*2*ATN2(SQRT((sin((radians(GeoLatitude - @lat2Compare)) / 2) * sin((radians(GeoLatitude - @lat2Compare)) / 2)) + (cos(radians(GeoLatitude)) * cos(radians(@lat2Compare)) * sin(radians(GeoLongitude - @long2Compare)/2) * sin(radians(GeoLongitude - @long2Compare)/2)))
    , SQRT(1-((sin((radians(GeoLatitude - @lat2Compare)) / 2) * sin((radians(GeoLatitude - @lat2Compare)) / 2)) + (cos(radians(GeoLatitude)) * cos(radians(@lat2Compare)) * sin(radians(GeoLongitude - @long2Compare)/2) * sin(radians(GeoLongitude - @long2Compare)/2)))
    ))) <= @radiusInKm

If you would like to perform the Haversine formula in C#,

    double resultDistance = 0.0;
    double avgRadiusOfEarth = 6371.392896; //Radius of the earth differ, I'm taking the average.

    //Haversine formula
    //distance = R * 2 * aTan2 ( square root of A, square root of 1 - A )
    //                   where A = sinus squared (difference in latitude / 2) + (cosine of latitude 1 * cosine of latitude 2 * sinus squared (difference in longitude / 2))
    //                   and R = the circumference of the earth

    double differenceInLat = DegreeToRadian(currentLatitude - latitudeToCompare);
    double differenceInLong = DegreeToRadian(currentLongitude - longtitudeToCompare);
    double aInnerFormula = Math.Cos(DegreeToRadian(currentLatitude)) * Math.Cos(DegreeToRadian(latitudeToCompare)) * Math.Sin(differenceInLong / 2) * Math.Sin(differenceInLong / 2);
    double aFormula = (Math.Sin((differenceInLat) / 2) * Math.Sin((differenceInLat) / 2)) + (aInnerFormula);
    resultDistance = avgRadiusOfEarth * 2 * Math.Atan2(Math.Sqrt(aFormula), Math.Sqrt(1 - aFormula));

DegreesToRadian is a function I custom created, its is a simple 1 liner of"Math.PI * angle / 180.0

My blog entry - SQL Haversine

Answer 3

Have a look at this.. has a javascript example as well.

Find Distance

Answer 4

Use the Great Circle Distance Formula.

Answer 5

If you are measuring distances less than (perhaps) 1 degree lat/long change, are looking for a very high performance approximation, and are willing to accept more inaccuracy than Haversine formula, consider these two alternatives:

(1) "Polar Coordinate Flat-Earth Formula" from Computing Distances:

a = pi/2 - lat1  
b = pi/2 - lat2  
c = sqrt( a^2 + b^2 - 2 * a * b * cos(lon2 - lon1) )   
d = R * c

(2) Pythagorean theorem adjusted for latitude, as seen in Ewan Todd's SO post:

d_ew = (long1 - long0) * cos(average(lat0, lat1))  
d_ns = (lat1 - lat0)  
d = sqrt(d_ew * d_ew + d_ns * d_ns)  

NOTES:
Compared to Ewan's post, I've substituted average(lat0, lat1) for lat0 inside of cos( lat0 ).

#2 is vague on whether values are degrees, radians, or kilometers; you will need some conversion code as well. See my complete code at bottom of this post.

#1 is designed to work well even near the poles, though if you are measuring a distance whose endpoints are on "opposite" sides of the pole (longitudes differ by more than 90 degrees?), Haversine is recommended instead, even for small distances.

I haven't thoroughly measured errors of these approaches, so you should take representative points for your application, and compare results to some high-quality library, to decide if the accuracies are acceptable. For distances less than a few kilometers my gut sense is that these are within 1% of correct measurement.

---

An alternative way to gain high performance (when applicable):

If you have a large set of static points, within one or two degrees of longitude/latitude, that you will then be calculating distances from a small number of dynamic (moving) points, consider converting your static points ONCE to the containing UTM zone (or to any other local Cartesian coordinate system), and then doing all your math in that Cartesian coordinate system.
Cartesian = flat earth = Pythagorean theorem applies, so distance = sqrt(dx^2 + dy^2).

Then the cost of accurately converting the few moving points to UTM is easily afforded.

---

CAVEAT for #1 (Polar): May be very wrong for distances less than 0.1 (?) meter. Even with double precision math, the following coordinates, whose true distance is about 0.005 meters, was given as "zero" by my implementation of Polar algorithm:

inputs:

    lon1Xdeg    16.6564465477996    double
    lat1Ydeg    57.7760262271983    double
    lon2Xdeg    16.6564466358281    double
    lat2Ydeg    57.776026248554 double

results:

Oblate spheroid formula:  
    0.00575254911118364 double
Haversine:
    0.00573422966122257 double
Polar:
    0

this was due to the two factors u and v exactly canceling each other:

    u   0.632619944868587   double
    v   -0.632619944868587  double

In another case, it gave a distance of 0.067129 m when the oblate spheroid answer was 0.002887 m. The problem was that cos(lon2 - lon1) was too close to 1, so cos function returned exactly 1.

Other than measuring sub-meter distances, the max errors (compared to an oblate spheroid formula) I found for the limited small-distance data I've fed in so far:

    maxHaversineErrorRatio  0.00350976281908381 double
    maxPolarErrorRatio  0.0510789996931342  double

where "1" would represent a 100% error in the answer; e.g. when it returned "0", that was an error of "1" (excluded from above "maxPolar"). So "0.01" would be an error of "1 part in 100" or 1%.

Comparing Polar error with Haversine error over distances less than 2000 meters to see how much worse this simpler formula is. So far, the worst I've seen is 51 parts per 1000 for Polar vs 4 parts per 1000 for Haversine. At about 58 degrees latitude.

---

Now implemented "Pythagorean with Latitude Adjustment".

It is MUCH more consistent than Polar for distances < 2000 m.
I originally thought the Polar problems were only when < 1 m,
but the result shown immediately below is quite troubling.

As distances approach zero, pythagorean/latitude approaches haversine. For example this measurement ~ 217 meters:

    lon1Xdeg    16.6531667510102    double
    lat1Ydeg    57.7751705615804    double
    lon2Xdeg    16.6564468739869    double
    lat2Ydeg    57.7760263007586    double

    oblate      217.201200413731
    haversine   216.518428601051
    polar       226.128616011973
    pythag-cos  216.518428631907
    havErrRatio 0.00314349925958048
    polErrRatio 0.041102054598393
    pycErrRatio 0.00314349911751603

Polar has a much worse error with these inputs; either there is some mistake in my code, or in Cos function I am running on, or I have to recommend not using Polar, even though most Polar measurements were much closer than this.

OTOH, Pythagorean, even with * cos(latitude) adjustment, has error that increases more rapidly than distance (ratio of max_error/distance increases for larger distances), so you need to carefully consider the maximum distance you will measure, and the acceptable error. In addition, it is not advisable to COMPARE two nearly-equal distances using Pythagorean, to decide which is shorter, as the error is different in different DIRECTIONS (evidence not shown).

Worst case measurements, errorRatio = Abs(error) / distance (Sweden; up to 2000 m):

    t_maxHaversineErrorRatio    0.00351012021578681 double
    t_maxPolarErrorRatio        66.0825360597085    double
    t_maxPythagoreanErrorRatio  0.00350976281416454 double

As mentioned before, the extreme polar errors are for sub-meter distances, where it could report zero instead of 6 cm, or report over 0.5 m for a distance of 1 cm (hence the "66 x" worst case shown in t_maxPolarErrorRatio), but there are also some poor results at larger distances. [Needs to be tested again with a Cosine function that is known to be highly accurate.]

Measurements taken in C# code in Xamarin.Android running on a Moto E4.

---

C# code:

    // x=longitude, y= latitude. oblate spheroid formula. TODO: From where?
    public static double calculateDistanceDD_AED( double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg )
    {
        double c_dblEarthRadius = 6378.135; // km
        double c_dblFlattening = 1.0 / 298.257223563; // WGS84 inverse
                                                      // flattening
        // Q: Why "-" for longitudes??
        double p1x = -degreesToRadians( lon1Xdeg );
        double p1y = degreesToRadians( lat1Ydeg );
        double p2x = -degreesToRadians( lon2Xdeg );
        double p2y = degreesToRadians( lat2Ydeg );

        double F = (p1y + p2y) / 2;
        double G = (p1y - p2y) / 2;
        double L = (p1x - p2x) / 2;

        double sing = Math.Sin( G );
        double cosl = Math.Cos( L );
        double cosf = Math.Cos( F );
        double sinl = Math.Sin( L );
        double sinf = Math.Sin( F );
        double cosg = Math.Cos( G );

        double S = sing * sing * cosl * cosl + cosf * cosf * sinl * sinl;
        double C = cosg * cosg * cosl * cosl + sinf * sinf * sinl * sinl;
        double W = Math.Atan2( Math.Sqrt( S ), Math.Sqrt( C ) );
        if (W == 0.0)
            return 0.0;

        double R = Math.Sqrt( (S * C) ) / W;
        double H1 = (3 * R - 1.0) / (2.0 * C);
        double H2 = (3 * R + 1.0) / (2.0 * S);
        double D = 2 * W * c_dblEarthRadius;

        // Apply flattening factor
        D = D * (1.0 + c_dblFlattening * H1 * sinf * sinf * cosg * cosg - c_dblFlattening * H2 * cosf * cosf * sing * sing);

        // Transform to meters
        D = D * 1000.0;

        // tmstest
        if (true)
        {
            // Compare Haversine.
            double haversine = HaversineApproxDistanceGeo( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg );
            double error = haversine - D;
            double absError = Math.Abs( error );
            double errorRatio = absError / D;
            if (errorRatio > t_maxHaversineErrorRatio)
            {
                if (errorRatio > t_maxHaversineErrorRatio * 1.1)
                    Helper.test();
                t_maxHaversineErrorRatio = errorRatio;
            }

            // Compare Polar Coordinate Flat Earth. 
            double polarDistanceGeo = ApproxDistanceGeo_Polar( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
            double error2 = polarDistanceGeo - D;
            double absError2 = Math.Abs( error2 );
            double errorRatio2 = absError2 / D;
            if (errorRatio2 > t_maxPolarErrorRatio)
            {
                if (polarDistanceGeo > 0)
                {
                    if (errorRatio2 > t_maxPolarErrorRatio * 1.1)
                        Helper.test();
                    t_maxPolarErrorRatio = errorRatio2;
                }
                else
                    Helper.dubious();
            }

            // Compare Pythagorean Theorem with Latitude Adjustment. 
            double pythagoreanDistanceGeo = ApproxDistanceGeo_PythagoreanCosLatitude( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
            double error3 = pythagoreanDistanceGeo - D;
            double absError3 = Math.Abs( error3 );
            double errorRatio3 = absError3 / D;
            if (errorRatio3 > t_maxPythagoreanErrorRatio)
            {
                if (D < 2000)
                {
                    if (errorRatio3 > t_maxPythagoreanErrorRatio * 1.05)
                        Helper.test();
                    t_maxPythagoreanErrorRatio = errorRatio3;
                }
            }
        }


        return D;
    }

    // As a fraction of the distance.
    private static double t_maxHaversineErrorRatio, t_maxPolarErrorRatio, t_maxPythagoreanErrorRatio;


    // Average of equatorial and polar radii (meters).
    public const double EarthAvgRadius = 6371000;
    public const double EarthAvgCircumference = EarthAvgRadius * 2 * PI;
    // CAUTION: This is an average of great circles; won't be the actual distance of any longitude or latitude degree.
    public const double EarthAvgMeterPerGreatCircleDegree = EarthAvgCircumference / 360;

    // Haversine formula (assumes Earth is sphere).
    // "deg" = degrees.
    // Perhaps based on Haversine Formula in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
    public static double HaversineApproxDistanceGeo(double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg)
    {
        double lon1 = degreesToRadians( lon1Xdeg );
        double lat1 = degreesToRadians( lat1Ydeg );
        double lon2 = degreesToRadians( lon2Xdeg );
        double lat2 = degreesToRadians( lat2Ydeg );

        double dlon = lon2 - lon1;
        double dlat = lat2 - lat1;
        double sinDLat2 = Sin( dlat / 2 );
        double sinDLon2 = Sin( dlon / 2 );
        double a = sinDLat2 * sinDLat2 + Cos( lat1 ) * Cos( lat2 ) * sinDLon2 * sinDLon2;
        double c = 2 * Atan2( Sqrt( a ), Sqrt( 1 - a ) );
        double d = EarthAvgRadius * c;
        return d;
    }

    // From https://stackoverflow.com/a/19772119/199364
    // Based on Polar Coordinate Flat Earth in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
    public static double ApproxDistanceGeo_Polar( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
    {
        double approxUnitDistSq = ApproxUnitDistSq_Polar(lon1deg, lat1deg, lon2deg, lat2deg, D);
        double c = Sqrt( approxUnitDistSq );
        return EarthAvgRadius * c;
    }

    // Might be useful to avoid taking Sqrt, when comparing to some threshold.
    // Threshold would have to be adjusted to match:  Power(threshold / EarthAvgRadius, 2)
    private static double ApproxUnitDistSq_Polar(double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
    {
        const double HalfPi = PI / 2; //1.5707963267949;

        double lon1 = degreesToRadians(lon1deg);
        double lat1 = degreesToRadians(lat1deg);
        double lon2 = degreesToRadians(lon2deg);
        double lat2 = degreesToRadians(lat2deg);

        double a = HalfPi - lat1;
        double b = HalfPi - lat2;
        double u = a * a + b * b;
        double dlon21 = lon2 - lon1;
        double cosDeltaLon = Cos( dlon21 );
        double v = -2 * a * b * cosDeltaLon;
        // TBD: Is "Abs" necessary?  That is, is "u + v" ever negative?
        //   (I think not; "v" looks like a secondary term. Though might be round-off issue near zero when a~=b.)
        double approxUnitDistSq = Abs(u + v);

        //if (approxUnitDistSq.nearlyEquals(0, 1E-16))
        //  Helper.dubious();
        //else if (D > 0)
        //{
        //  double dba = b - a;
        //  double unitD = D / EarthAvgRadius;
        //  double unitDSq = unitD * unitD;
        //  if (approxUnitDistSq > 2 * unitDSq)
        //      Helper.dubious();
        //  else if (approxUnitDistSq * 2 < unitDSq)
        //      Helper.dubious();
        //}

        return approxUnitDistSq;
    }

    // Pythagorean Theorem with Latitude Adjustment - from Ewan Todd - https://stackoverflow.com/a/1664836/199364
    // Refined by ToolmakerSteve - https://stackoverflow.com/a/53468745/199364
    public static double ApproxDistanceGeo_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
    {
        double approxDegreesSq = ApproxDegreesSq_PythagoreanCosLatitude( lon1deg, lat1deg, lon2deg, lat2deg );
        // approximate degrees on the great circle between the points.
        double d_degrees = Sqrt( approxDegreesSq );
        return d_degrees * EarthAvgMeterPerGreatCircleDegree;
    }

    public static double ApproxDegreesSq_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg )
    {
        double avgLatDeg = average( lat1deg , lat2deg );
        double avgLat = degreesToRadians( avgLatDeg );

        double d_ew = (lon2deg - lon1deg) * Cos( avgLat );
        double d_ns = (lat2deg - lat1deg);
        double approxDegreesSq = d_ew * d_ew + d_ns * d_ns;
        return approxDegreesSq;
    }

Answer 6

here is a fiddle with finding locations / near locations to long/lat by given IP:

http://jsfiddle.net/bassta/zrgd9qc3/2/

And here is the function I use to calculate the distance in straight line:

function distance(lat1, lng1, lat2, lng2) {
        var radlat1 = Math.PI * lat1 / 180;
        var radlat2 = Math.PI * lat2 / 180;
        var radlon1 = Math.PI * lng1 / 180;
        var radlon2 = Math.PI * lng2 / 180;
        var theta = lng1 - lng2;
        var radtheta = Math.PI * theta / 180;
        var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
        dist = Math.acos(dist);
        dist = dist * 180 / Math.PI;
        dist = dist * 60 * 1.1515;

        //Get in in kilometers
        dist = dist * 1.609344;

        return dist;
    }

It returns the distance in Kilometers