Smoky Cogs

Don’t you think it would be really great to be able to work out the precise bearing of a place in relation to where you are now? The calculation for this is rather simple if you have the geographical coordinates of the two points that you would like to find the bearing of.

The code below calculates the angle between the two coordinates. The last function call to ConvertToBearing puts the value into the right range. The atan function returns a value back (converted to degrees), which is in the range of -180 to 180 degrees, whereas, bearing is normally given in the range 0 – 360 degrees.

namespace Geodesy
{
    public class Geodesy
    {
        public static double RadToDeg(double radians)
        {
            return radians * (180 / Math.PI);
        }

        public static double DegToRad(double degrees)
        {
            return degrees * (Math.PI / 180);
        }

        public static double Bearing(double lat1, double long1, double lat2, double long2)
        {
            //Convert input values to radians
            lat1 = Geodesy.DegToRad(lat1);
            long1 = Geodesy.DegToRad(long1);
            lat2 = Geodesy.DegToRad(lat2);
            long2 = Geodesy.DegToRad(long2);

            double deltaLong = long2 - long1;

            double y = Math.Sin(deltaLong) * Math.Cos(lat2);
            double x = Math.Cos(lat1) * Math.Sin(lat2) -
                    Math.Sin(lat1) * Math.Cos(lat2) * Math.Cos(deltaLong);
            double bearing = Math.Atan2(y, x);
            return Geodesy.ConvertToBearing(Geodesy.RadToDeg(bearing));
        }

        public static double ConvertToBearing(double deg)
        {
            return (deg + 360) % 360;
        }
    }
}

Previously I showed how to find the distance between coordinates and the bearing, but what if you want to find the midpoint of two coordinates.

The code in C# is rather straightforward – that is, if you know a little bit of trigonometry.

namespace Geodesy
{
    public class Geodesy
    {
        public static double RadToDeg(double radians)
        {
            return radians * (180 / Math.PI);
        }

        public static double DegToRad(double degrees)
        {
            return degrees * (Math.PI / 180);
        }

        public static void MidpointCoords(double lat1, double long1, double lat2, double long2, ref double midpointLat, ref double midpointLong)
        {
            //Convert input values to radians
            lat1 = Geodesy.DegToRad(lat1);
            long1 = Geodesy.DegToRad(long1);
            lat2 = Geodesy.DegToRad(lat2);
            long2 = Geodesy.DegToRad(long2);
            double deltaLat = lat2 - lat1;
            double deltaLong = long2 - long1;

            double Bx = Math.Cos(lat2) * Math.Cos(deltaLong);
            double By = Math.Cos(lat2) * Math.Sin(deltaLong);
            midpointLat = Geodesy.RadToDeg(Math.Atan2(Math.Sin(lat1) + Math.Sin(lat2),
               Math.Sqrt((Math.Cos(lat1) + Bx) * (Math.Cos(lat1) + Bx) + By * By)));
            midpointLong = Geodesy.RadToDeg(long1 + Math.Atan2(By, Math.Cos(lat1) + Bx));

        }

    }
}

With the advent of GPS systems, adn mapping systems like Google Earth, it has become rather useful to be able to calculate the distance (as the crow flies) between two points on the earth.

One method to do so is known as the Haversine formula.

Taking the longitude and latitude of two points, and then converting them to radians, since the C# trig functions take values in radians, not in degrees.

We then take the difference of latitude and longitude of the two points.

Applying the ‘Haversine’ formula needed and then multiplying by the radius of the earth gives us the distance in kilometres of the two points.
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
Distance = (radius of earth) * c

This method of calculating the distance is not entirely accurate though, but will suffice for all but the most exacting of applications. For every kilometer, there is an average error of 3 meters, or roughly a 0.3-0.5% error. This error creeps in because the formula assumes a perfectly spherical earth, whereas the earth is slightly ellipsoid, so the radius around the equator is slightly bigger than the radius around the poles.

namespace Geodesy
{
    public class Geodesy
    {
        public static double RadToDeg(double radians)
        {
            return radians * (180 / Math.PI);
        }

        public static double DegToRad(double degrees)
        {
            return degrees * (Math.PI / 180);
        }


        //Calculate distance between two coordinates useing the Haversine formula
        public static double DistanceBetweenCoords(double lat1, double long1, double lat2, double long2)
        {
            //Convert input values to radians
            lat1 = Geodesy.DegToRad(lat1);
            long1 = Geodesy.DegToRad(long1);
            lat2 = Geodesy.DegToRad(lat2);
            long2 = Geodesy.DegToRad(long2);

            double earthRadius = 6371; // km
            double deltaLat = lat2 - lat1;
            double deltaLong = long2 - long1;
            double a = Math.Sin(deltaLat / 2.0) * Math.Sin(deltaLat / 2.0) +
                    Math.Cos(lat1) * Math.Cos(lat2) *
                    Math.Sin(deltaLong / 2.0) * Math.Sin(deltaLong / 2.0);
            double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
            double distance = earthRadius * c;

            return distance;
        }

}