Geographic Distance and Azimuth Calculations

Introduction

Many applications, such as navigation and radio frequency engineering, require a thorough understanding of geographic calculations. Some very natural questions seem to come up in a variety of disciplines: How far apart are two points on the Earth? What direction do I need to go to reach a particular point? If I go in a particular direction for a certain distance, where will I end up?

Visualization makes these calculations immensely easier, and to visualize you need to come up with an accurate model. As it turns out, there are two common approaches for modelling the surface of the Earth: spherical and ellipsoidal. Another, more accurate, model is called the geoid. The geoid is a complex surface where each point on the surface has the same gravitational pull. The shape of the geoid does not lend itself well to geometric calculations (and new research and measurements are constantly refining the geoid), so people generally stick with either the spherical model or the ellipsoid model. The spherical model can be very accurate under certain stringent conditions; however, the ellipsoid model is generally a very accurate model everywhere. You can think of either model as the mean sea level. So elevations, such as those on a contour map, are generally given as height above the ellipsoid. Both spherical and ellipsoid models have symmetry that allow you to do calculations, but that symmetry also means that people have to agree on a common starting point for the model. The starting "reference" point is called a datum and there are many different datums. Transforming between datums can be very complicated depending on how you do it, and those transformations are just outside the scope of this article (maybe another article will cover datums). So, the rest of this article assumes that we are working within some particular datum and there is no need to transfer the coordinates in this datum into coordinates of another datum.

The good news is that we can solve a lot of geographic problems in the spherical model with a few simple mathematical tools. Another important aspect of the spherical model is that, in terms of visualization, it covers just about everything we need; the ellipsoid model can be visualized as a refinement of the spherical model. This approach works really well because, in terms of percentages, the Earth is very close to a sphere.

Overview

This article describes each model in some depth and provides solutions to the following common geographic problems:

  1. Spherical Model
    1. Calculate path length given starting and ending coordinates (a pair of latitude/longitude).
    2. Calculate path direction (azimuth) given starting and ending coordinates.
    3. Calculate end point (latitude/longitude) given a starting point, distance, and azimuth.
    4. Calculate the azimuth and elevation angle to a geostationary satellite (where to point a satellite dish).
    5. Calculate the intersection of two paths given starting and ending coordinates for each path.
  2. Ellipsoid Model
    1. Calculate path length along a meridian given starting and ending coordinates.
    2. Calculate azimuth and path length (the geodetic forward and inverse problems).

Spherical Model

The spherical model is simple in mathematical terms because of its symmetry: every point on the surface is equidistant from the center, it's difficult to imagine more symmetry. This fact has a number of very helpful consequences that can be summed up in the following statment: Geodesic paths between two points on a sphere are great circles. A geodesic path is simply the shortest path between two points along the surface. Of course, it would be shorter to go straight through the Earth between the two points, but that is generally not possible for us surface dwellers. A great circle is just like every other circle with the additional contraint that its center lies at the center of the sphere. The following table summarizes some of the mathematical tools that are available for analyzing the spherical model:

Law of Cosines (spherical) cos(b) = sin(a)*sin(c) + cos(a)*cos(c)*cos(B)
Law of Sines (spherical) sin (B)/sin(b) = sin(A)/sin(a) = sin(C)/sin(c)
Circular Arc Length length = (radius) * (angular separation) = r * Q
Pythagorean Theorem (3D) d = sqrt(x*x + y*y + z*z)
Law of Cosines (plane) c^2 = a^2 + b^2 - 2*a*b*cosg
Table 1: Mathematical tools available for the spherical model.

As with all mathematical formulas, it's important to understand where these formulas apply. For instance, lines of latitude are NOT great circles (except for the equator, which is a line of latitude and a great circle) and so laws of Cosines and Sines for spherical triangles do not apply to lines of latitude. Lines of latitude are in fact circles with their centers lying along the polar axis of the Earth, not necessarily at the center of the Earth. Lines of longitude are great circles. Figure 1 depicts a generic great circle path.


Figure 1: A great circle path (arc b is the path and the angular separation between the end points).

For great circle paths, the start and end points as well as all the points along the path and the center of the Earth lie in the same plane. Use Figure 2 to visualize the geometry of a great circle path/plane.

Figure 2: A rectangular plane intersecting a great circle path and the center of the Earth (arc b is the path and the angular separation of the end points).

Problem 1A. Calculate path length given starting and ending coordinates

Calculating distance for a spherical model is very straightforward using the Law of Cosines and the Law of Sines for spherical triangles. You may remember the Law of Cosines and the Law of Sines for planar triangles; well, it's slightly different for spherical triangles. In Figure 2, the triangle with sides {a, b, c} and angles {A, B, C} is defined by two end points and the North Pole. Given the latitude and longitude for the end points, we want to solve for side b and angle A, so we start with the Law of Cosines for spherical triangles:

cos (b) = cos (a) * cos (c) + sin (a) * sin (c) * cos (B)

where B = lon2 - lon1, and c = 90 - lat1, and a = 90 - lat2, and substituting these above leads to

cos (b) = cos (90 - lat2) * cos (90 - lat1) + sin (90 - lat2) * sin (90 - lat1) * cos (lon2 - lon1), and solving for b

b = arccos ( cos (90 - lat2) * cos (90 - lat1) + sin (90 - lat2) * sin (90 - lat1) * cos (lon2 - lon1) )

Because b equals the arccos of something, b is an angle. It's actually the angular distance between the two points (see Figure 3), and for circles the arc length is:

arc length =( radius )* ( angular distance (in radians)), so finally

distance = ( Earth Radius ) * arccos ( cos (90 - lat2) * cos (90 - lat1) + sin (90 - lat2) * sin (90 - lat1) * cos (lon2 - lon1) )

If you use miles for the Earth Radius, you will get the distance in miles for a result; using kilometers yields the answer in kilometers, and so on.

One (1st edition) reader provided an alternate distance formula that is better for small angles. I'm not sure who first proposed using this formula, but it seems to work well.

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

If there were a need to calculate the distance between two points that are very close together, it would be tempting to use the brute force approach of the Pythagorean Theorem. Essentially, it means calculating the {x,y,z} position of each point (in whatever units suit you) and using distance = sqrt(x*x + y*y +z*z). The beauty of this formula is it has no undefined arguments. Also, when you are dealing with small distances, the elevation difference between the end points may make a bigger difference than the curvature of the Earth. For instance, if you have the geographic positions and elevations of the terminals of a very long gondola (maybe from a GPS receiver) and you want to calculate how far it has to travel (and for some reason it's impossible to measure the length of cable used). You could use the following formulas to calculate {x,y,z} and plug them into the 3D distance formula:

x1 = (R + h1) * [ cos(lat1) * cos(lon1) ] ; y1 = (R + h1) * [ sin(lat1)]; z1 = (R + h1) * [ cos(lat1) * sin(lon1) ]

x2 = (R + h2) * [ cos(lat2) * cos(lon2) ] ; y2 = (R + h2) * [ sin(lat2)]; z2 = (R + h2) * [ cos(lat2) * sin(lon2) ]

dx = (x2 - x1); dy = (y2 - y1); dz = (z2 - z1)

distance = sqrt(dx*dx + dy*dy + dz*dz) where h1 and h2 are the elevation (above the mean sea level) of the two end points and R is the radius of the Earth.

This approach would be very accurate for distances less than 10 miles. The inaccuracy of this approach, obviously, is that it assumes that you can travel in a perfectly straight line between the points. For longer distances (actually for all paths to some degree), the curvature of the Earth would prevent this from happening.

Problem 1B. Calculate path direction (azimuth) given starting and ending coordinates

Using the results of Problem 1A and the Law of Sines, find the azimuth from (lat1,lon1) to (lat2,lon2). By the Law of Sines:

sin (A) / sin (a) = sin (B) / sin (b) = sin (C) / sin (c). We want to calculate A and we just found b in Problem 1A, so we have everything we need:

sin (A) = sin (a) * sin (B) / sin (b); or A = arcsin ( sin (90 - lat2) * sin (lon2 - lon1) / sin (b) )

Once you find A, the azimuth can be determined based on the relationship of the end points. For Figure 2, A is equal to the azimuth. The sample project provides a method for determining the azimuth from A depending on the relationship of the end points.

Problem 1C. Calculate end point (latitude/longitude) given a starting point, distance, and azimuth

Given {lat1, lon1, distance, azimuth} calculate {lat2, lon2}. First, work backwards (relative to Problem 1A) and find b from the distance by dividing by the Earth radius. b = distance / (Earth Radius) making sure distance and (Earth Radius) are the same units so that we end up with b in radians. Knowing b, calculate a using a = arccos(cos(b)*cos(90 - lat1) + sin(90 - lat1)*sin(b)*cos(azimuth))—basically taking the arc cosine of the Law of Cosines for a. From a, we can get lat2, so the only item remaining is to figure lon2; we can get that if we know B. Calculate B using B = arcsin(sin(b)*sin(azimuth)/sin(a)). Then finally, lat2 = 90 - a and lon2 = B + lon1. Essentially, we just worked Problem 1A backwards.

Problem 1D. Calculate the azimuth and elevation angle to a geostationary satellite (where to point a satellite dish)

Many satellite-direct TV services use geostationary satellites to transmit their signals to peoples' homes. Geostationary satellties move in nearly circular orbits approximately above the equator. Because a geostationary satellite keeps pace with the rotation of the Earth, it appears fixed in space to anyone on the surface of the Earth (anyone who is not moving). So, the geometry of the orbit and the fixed nature of the users enables home users to point their antennas (small satellite dish) in a fixed direction. But, what is the direction to the satellite and how much tilt above the horizon should a dish have? Believe it or not, these are simple geographic calculations! All of the geostationary satellites lie very nearly above the equator, so we know lat2. The question is: What is the longitude? Well, it varies depending on which satellite you want to reach. The service providers want to reach as many households as possible, so for the western hemisphere they want to position the satellite above the northern part of South America. From there, customers in most of North and South America can see the satellite and pay them for the service. The list below shows some of the longitudes of the major satellite TV stations:

DirectTV 1 109.8 W
DirectTV 1R 101.2 W
DirectTV 2 100.8 W
DirectTV 3 100.8 W
DirectTV 4S 101.2 W
DirectTV 6 118.8 W
EchoStar 5/Dish Network 110.0 W
EchoStar 2 119.0 W
EchoStar 4 119.5 W
EchoStar 6 119.0 W
EchoStar 7 119.0 W
Table 2: Direct broadcast service satellite locations.

The altitude of geostationary satellites is about 35,800 km or 22,300 miles.

Figure 3: Geostationary satellite orbit (over the Atlantic Ocean) and pointing geometry

Figure 4 shows a more detailed visualization of the path from the satellite to your receiver. For this type of problem, we will use a combination of spherical geometry and plane geometry. The plane geometry will be used to analyze the shaded triangle in Figure 4.

Figure 4: Detailed path geometry from a geostationary satellite to an Earth-based receiver (r = receiver for this figure, R = Earth Radius).

Okay, let's get to it. So, we have a receiver r point at {lat1, lon1} and a transmitter at { lat2 = 0.0 (equator), lon2} and we want to find the azimuth, and elevation (tilt) angle. The azimuth is precisely the same that we computed in Problem 1B with lat2 = 0.0 and lon2 = longitude of the geostationary satellite. I'll repeat those calculations here so you don't have to look up too much:

b = acos(cos(lon2-lon1)*cos(lat1)) and A = arcsin ( sin (90 - lat2) * sin (lon2 - lon1) / sin (b) ), and from A you find the azimuth in the same way described above (in Problem 1B). To find the elevation angle, we need to take a closer look at the triangle in Figure 4c. The angle of the triangle at point r is greater than 90 degrees (if the receiver can see the satellite). If the angle at r were 90 degrees, the satellite would be right on the horizon, so we will measure the tilt angle as the amount of angle r above 90 degrees. Using the Law of Cosines for Plane triangles, we get: (R+h)^2 = d^2 + R^2 - 2*R*d*cos(r), so r = arccos((d^2 + R^2 - (R+h)^2)/(2*R*d)), where R is the radius of the Earth, h the altitude of the geostationary satellite, and d the distance from the receiver to the satellite. Simplifying and putting things in terms of values that we know:

d = sqrt(R*R + (R+h)*(R+h) - 2*R*(R+h)*cos(b)) and elevation angle (tilt up from horizon) = arccos((R+h)*sin(b)/d).

Problem 1E. Calculate the intersection of two paths given the starting and ending coordinates for each path

Unfortunately, calculating the intersection of two great circle paths requires slightly more mathematical tools than what appear in Table 1. My goal in writing this article is to keep the math to a bare minimum and justify as much as possible with visualization rather than formulas; however, the easiest way (I know of) to calculate the intersection of two paths makes use of some basic vector analysis. For those who need it, the next paragraph provides a brief vector review.

First, a very brief review of vector arithmetic: A vector is a combined representation of a length and direction. For example, the vectors v1 = {1, 1, 1} and v2 = {2, 2, 2} have the same direction (because the ratio of the components are the same); however, v2 is longer (because the values are larger). The length of a vector is equal to the square root of the sum of its components-squared, for example: Length(v1) = sqrt(v1[0]*v1[0] + v1[1]*v[1] + v1[2]*v1[2]). Vectors can be multiplied in two distinct ways: the scalar (dot) product or the vector (cross) product. The great thing about the vector product is that this form of multiplication produces another vector that is guaranteed to be perpendicular (okay, "normal" if you want to be picky) to the original two vectors. This fact can really help since Great Circle paths lie in a plane.

If vectors are too confusing, just keep in mind that we have a simple formula that takes two vectors as inputs and produces another vector that is guaranteed perpendicular to both of the input vectors. Also keep in mind that it takes three distinct points to define a plane.


Figure 5: Geometry of intersection (a) Earth and both planes and normals, (b) plane and normal of path A, (c) plane and normal of path B, (d) both planes, both normals and the intersection vector.

Let's take the paths one at a time. Path A depicted in Figure 5b: the two points on the plane are the end points of the path. If we take the center of the earth, and the starting and ending points, we have a plane. To define the plane, we need the normal vector, also known as the normal. We can easily get that by taking a vector from the center of the Earth to the starting point for path A (vector1), and a vector from the center of the Earth to the ending point of path A (vector2). Both vector1 and vector2 are "in the plane." So, if we take the vector product of vector1 and vector2, we get the normal for path A (normalA). Figure 5c is the same picture for path B. Repeating the procedure for path B, we get the normal for path B (normalB).

So, normalA is perpendicular to every vector in the plane of path A, and normalB is perpendicular to every vector in the plane of path B. The only vector that the two planes share is the intersection vector. Because the intersection is a vector in the plane of path A, and also in the plane of path B, it must be perpendicular to both normals. So, if we take the vector product of the normals, we are sure to get the intersection vector! Wow, I just realized how difficult it is to put that kind of analysis into written words. For those who prefer to read code, the following snippet calculates the point where the intersection vector passes through the Earth's surface; in other words, a pair of latitude and longitude coordinates on opposite sides of the Earth.

For the intersect problem, it gets a little tricky to keep track of all the coordinates, so let me explain. { lat1A, lon1A, lat1B, lon1B } are the end-points of path 1. { lat2A, lon2A, lat2B, lon2B } are the end-points of path 2. { lat3A, lon3A, lat3B, lon3B } are the points where the intersection vector passes through the surface of the Earth. Keep in mind the the planes of the paths will always intersect, but the paths themselves may not intersect. In case of path-intersection, I put the path intersection point into { lat3A, lon3A } and the intersection coordinate on the opposite side of the Earth into { lat3B, lon3B }. For the case of non-intersecting paths, { lat3A, lon3A } and { lat3B, lon3B } are just points where the intersection vector of the planes touch the surface of the Earth, and neither are inside the path segments. The return value (bool) is true if the segments intersect and false if they do not intersect, either way the plane intersection points are calculated.

namespace GEO {
   const double PI                = 3.14159265359;
   const double TWOPI             = 6.28318530718;
   const double DE2RA             = 0.01745329252;
   const double RA2DE             = 57.2957795129;
   const double ERAD              = 6378.135;
   const double ERADM             = 6378135.0;
   const double AVG_ERAD          = 6371.0;
   const double FLATTENING        = 1.0/298.257223563;
                                  // Earth flattening (WGS '84)
   const double EPS               = 0.000000000005;
   const double KM2MI             = 0.621371;
   const double GEOSTATIONARY_ALT = 35786.0;    // km
}

bool GCIntersectSegment(double lat1A, double lon1A, double lat1B,
                        double lon1B, double lat2A, double lon2A,
                        double lat2B, double lon2B, double& lat3A,
                        double& lon3A, double& lat3B, double& lon3B)
{
   bool isInside = false;

   double v1[3], v2[3], v3a[3], v3b[3], n1[3], n2[3];
   double m;

   double d1 = ApproxDistance(lat1A, lon1A, lat1B, lon1B);
   double d2 = ApproxDistance(lat2A, lon2A, lat2B, lon2B);

   //
   // for path 1, setting up my 2 vectors, v1 is vector
   // from center of the Earth to point A, v2 is vector
   // from center of the Earth to point B.
   //
   v1[0] = cos(lat1A * GEO::DE2RA) * cos(lon1A * GEO::DE2RA);
   v1[1] = sin(lat1A * GEO::DE2RA);
   v1[2] = cos(lat1A * GEO::DE2RA) * sin(lon1A * GEO::DE2RA);

   v2[0] = cos(lat1B * GEO::DE2RA) * cos(lon1B * GEO::DE2RA);
   v2[1] = sin(lat1B * GEO::DE2RA);
   v2[2] = cos(lat1B * GEO::DE2RA) * sin(lon1B * GEO::DE2RA);

   //
   // n1 is the normal to the plane formed by the three points:
   // center of the Earth, point 1A, and point 1B
   //
   n1[0] = (v1[1]*v2[2]) - (v1[2]*v2[1]);
   n1[1] = (v1[2]*v2[0]) - (v1[0]*v2[2]);
   n1[2] = (v1[0]*v2[1]) - (v1[1]*v2[0]);

   //
   // for path 2, setting up my 2 vectors, v1 is vector
   // from center of the Earth to point A, v2 is vector
   // from center of the Earth to point B.
   //
   v1[0] = cos(lat2A * GEO::DE2RA) * cos(lon2A * GEO::DE2RA);
   v1[1] = sin(lat2A * GEO::DE2RA);
   v1[2] = cos(lat2A * GEO::DE2RA) * sin(lon2A * GEO::DE2RA);

   v2[0] = cos(lat2B * GEO::DE2RA) * cos(lon2B * GEO::DE2RA);
   v2[1] = sin(lat2B * GEO::DE2RA);
   v2[2] = cos(lat2B * GEO::DE2RA) * sin(lon2B * GEO::DE2RA);

   //
   // n1 is the normal to the plane formed by the three points:
   // center of the Earth, point 2A, and point 2B
   //
   n2[0] = (v1[1]*v2[2]) - (v1[2]*v2[1]);
   n2[1] = (v1[2]*v2[0]) - (v1[0]*v2[2]);
   n2[2] = (v1[0]*v2[1]) - (v1[1]*v2[0]);

   //
   // v3 is perpendicular to both normal 1 and normal 2, so
   // it lies in both planes, so it must be the line of
   // intersection of the planes. The question is: does it
   // go towards the correct intersection point or away from
   // it.
   //
   v3a[0] = (n1[1]*n2[2]) - (n1[2]*n2[1]);
   v3a[1] = (n1[2]*n2[0]) - (n1[0]*n2[2]);
   v3a[2] = (n1[0]*n2[1]) - (n1[1]*n2[0]);

   //
   // want to make v3a a unit vector, so first have to get
   // magnitude, then each component by magnitude
   //
   m = sqrt(v3a[0]*v3a[0] + v3a[1]*v3a[1] + v3a[2]*v3a[2]);
   v3a[0] /= m;
   v3a[1] /= m;
   v3a[2] /= m;

   //
   // calculating intersection points 3A & 3B. A & B are
   // arbitrary designations right now, later we make A
   // the one close to, or within, the paths.
   //
   lat3A = asin(v3a[1]);
   if ((lat3A > GEO::EPS) || (-lat3A > GEO::EPS))
      lon3A = asin(v3a[2]/cos(lat3A));
   else
      lon3A = 0.0;

   v3b[0] = (n2[1]*n1[2]) - (n2[2]*n1[1]);
   v3b[1] = (n2[2]*n1[0]) - (n2[0]*n1[2]);
   v3b[2] = (n2[0]*n1[1]) - (n2[1]*n1[0]);

   m = sqrt(v3b[0]*v3b[0] + v3b[1]*v3b[1] + v3b[2]*v3b[2]);
   v3b[0] /= m;
   v3b[1] /= m;
   v3b[2] /= m;

   lat3B = asin(v3b[1]);
   if ((lat3B > GEO::EPS) || (-lat3B > GEO::EPS))
      lon3B = asin(v3b[2]/cos(lat3B));
   else
      lon3B = 0.0;

   //
   // get the distances from the path endpoints to the two
   // intersection points. these values will be used to determine
   // which intersection point lies on the paths, or which one
   // is closer.
   //
   double d1a3a = ApproxDistance(lat1A, lon1A, lat3A, lon3A);
   double d1b3a = ApproxDistance(lat1B, lon1B, lat3A, lon3A);
   double d2a3a = ApproxDistance(lat2A, lon2A, lat3A, lon3A);
   double d2b3a = ApproxDistance(lat2B, lon2B, lat3A, lon3A);

   double d1a3b = ApproxDistance(lat1A, lon1A, lat3B, lon3B);
   double d1b3b = ApproxDistance(lat1B, lon1B, lat3B, lon3B);
   double d2a3b = ApproxDistance(lat2A, lon2A, lat3B, lon3B);
   double d2b3b = ApproxDistance(lat2B, lon2B, lat3B, lon3B);

   if ((d1a3a < d1) && (d1b3a < d1) && (d2a3a < d2)
                    && (d2b3a < d2))
   {
      isInside = true;
   }
   else if ((d1a3b < d1) && (d1b3b < d1) && (d2a3b < d2)
                         && (d2b3b < d2))
   {
      //
      // 3b is inside the two paths, so swap 3a & 3b
      //
      isInside = true;

      m = lat3A;
      lat3A = lat3B;
      lat3B = m;
      m = lon3A;
      lon3A = lon3B;
      lon3B = m;
   }
   else
   {
      //
      // figure out which one is closer to the path
      //
      d1 = d1a3a + d1b3a + d2a3a + d2b3a;
      d2 = d1a3b + d1b3b + d2a3b + d2b3b;

      if (d1 > d2)
      {
         //
         // Okay, we are here because 3b {lat3B,lon3B} is closer to
         // the paths, so we need to swap 3a & 3b. the other case is
         // already the way 3a & 3b are organized, so no need to swap
         //
         m = lat3A;
         lat3A = lat3B;
         lat3B = m;
         m = lon3A;
         lon3A = lon3B;
         lon3B = m;
      }
   }

   //
   // convert the intersection points to degrees
   //
   lat3A *= GEO::RA2DE;
   lon3A *= GEO::RA2DE;
   lat3B *= GEO::RA2DE;
   lon3B *= GEO::RA2DE;

   return isInside;
}

Additional Code for Great Circle Calculations

double GCDistance(double lat1, double lon1, double lat2,
                  double lon2)
{
   lat1 *= GEO::DE2RA;
   lon1 *= GEO::DE2RA;
   lat2 *= GEO::DE2RA;
   lon2 *= GEO::DE2RA;
   double d = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1 -
                                                            lon2);
   return (GEO::AVG_ERAD * acos(d));
}

double GCAzimuth(double lat1, double lon1, double lat2, double lon2)
{
   double result = 0.0;

   INT32 ilat1 = (INT32)(0.50 + lat1 * 360000.0);
   INT32 ilat2 = (INT32)(0.50 + lat2 * 360000.0);
   INT32 ilon1 = (INT32)(0.50 + lon1 * 360000.0);
   INT32 ilon2 = (INT32)(0.50 + lon2 * 360000.0);

   lat1 *= GEO::DE2RA;
   lon1 *= GEO::DE2RA;
   lat2 *= GEO::DE2RA;
   lon2 *= GEO::DE2RA;

   if ((ilat1 == ilat2) && (ilon1 == ilon2))
   {
      return result;
   }
   else if (ilon1 == ilon2)
   {
      if (ilat1 > ilat2)
         result = 180.0;
   }
   else
   {
      double c = acos(sin(lat2)*sin(lat1) +
                 cos(lat2)*cos(lat1)*cos((lon2-lon1)));
      double A = asin(cos(lat2)*sin((lon2-lon1))/sin(c));
      result = (A * GEO::RA2DE);

      if ((ilat2 > ilat1) && (ilon2 > ilon1))
      {
      }
      else if ((ilat2 < ilat1) && (ilon2 < ilon1))
      {
         result = 180.0 - result;
      }
      else if ((ilat2 < ilat1) && (ilon2 > ilon1))
      {
         result = 180.0 - result;
      }
      else if ((ilat2 > ilat1) && (ilon2 < ilon1))
      {
         result += 360.0;
      }
   }

   return result;
}

As you can see, most of the code for calculating azimuths has to do with avoiding undefined arccos() situations and transforming A into Azimuth. The real calculations are only three lines long.

For all these problems, keep in mind that the standard C math library expects radians for arguments to sin() and cos(), and that acos() and asin() produce results in radians whereas latitude, longitude, and azimuth are usually expressed in degrees.

Ellipsoid Model

The geometry of the ellipsoid model is relatively simple to visualize by imagining cross sections of the Earth. A horizontal cross section (cutting through the Earth at the equator) would produce a circle of radius 6,378,137 meters. As the cross sections become more inclined with the equator, they become more elliptical (less circular). A vertical cross section that passes through the poles and the equator would be an ellipse with one axis (semi-major) radius equal to the equatorial radius (same as before) and the other axis (semi-minor, passing through the poles) with radius approximately 6,356,752 meters. This geometry can be conveniently described by two numbers: the equatorial radus ( a ) and the flattening ( f ). To calculate the polar radius, you simply multiply the equatorial radius times one minus the flattening. For example:

polar radius = equatorial radius * ( 1 - f ) or b = a * ( 1 - f )

Since the early 1800s, surveyors and geodisists have tried to estimate the flattening. Table 3 shows some of the common ellipsoids. Modern satellite-based measurements are far more accurate than previous estimates.

Ellipsoid Semi-major axis in meters Inverse flattening (1 / f )
Airy 1830 6,377,563.396 299.324964
Clarke 1866 6,378,206.4 294.978698
World Geodetic Survey (WGS72) 6,378,135 298.26
Geodetic Reference System 1980 (GRS80) 6,378,137 298.257222101
World Geodetic Survey (WGS84) 6,378,137 298.257223563
Table 3: Common Earth Ellipsoids

The technical name for a minimum-distance path is geodesic. Some geodesics are natural and intuitive while others can be quite complex. For example, most third graders know that the minimum distance path on a flat plane is a straight line. On the surface of a sphere, geodesics lie along great circles (i.e. a circle with its center at the center of the sphere). An ellipsoid geodesic is a minimum-distance path on an ellipsoid, and you might be tempted to think it's a great-ellipse (an ellipse where the center of the ellipse is located at the center of the Earth). If you think it's a great-ellipse, you are right in some cases; in other cases, it is not an ellipse but an ill-defined oval of some sort. This is unfortunate because even though solving for the arc length of an elliptical segment is painful, it's far easier than solving for (and visualizing) the differential geometry involved in ellipsoid geodesics. That's right, differential geometry...so much for our simple mathematical toolbox...it was fun while it lasted.

Problem 2A.: Calculate path length along a meridian given starting and ending coordinates

Try to visualize a path that runs along the equator, so lat1 = 0.0, and lat2 = 0.0. For the sake of simplicity, visualize a relatively short path, say 2 degrees arc length. The arc is a circlar-arc and calculating the distance and azimuth are trivial exercises using the toolbox from the spherical model (azimuth is simple zero, of course). Now, take point 2 and move it north of point 1 until they lie along the same meridian (longitude) except now lat2 > 0.0, say 2 degrees. In this case, the arc lies along an ellipse that runs through the North and South Poles. The semi-major axis of the ellipse is the equatorial radius ( a )of whatever ellipsoid we want to use, and the semi-minor axis is the polar radius, or as stated above b = a * (1 - f). This is the case we are dealing with for Problem 2A. In both cases, an equatorial arc and a polar arc, the center of the ellipse (remember a cicle is also an ellipse) is located at the center of the Earth. Now, just to complete the visualization, imagine point 1 and 2 at the same latitude again, this time let's say 60 degrees and separated by 2 degrees of longitude. In this case, the geodesic arc connecting the two points still lies in a plane, but unfortunately that plane does not contain the center of the Earth. This is the break-down in symmetry that makes it a really tough problem (Problem 2B). Okay, more on that in 2B; back to the case of the arc that lies along a meridian.

So, for this case we are dealing with an ellipse, an elliptical segment to be precise. Unfortunately, there is no handy dandy formula for the arc length of an ellipse, though many have tried and in the process come up with some terrific approximations. The basic approach for calculating the segment length of an ellipse is to break it down into tiny little segments, so tiny that each segment can be treated as a straight line. Then, take those straight lines and use them to form right triangles where the tiny segment is the hypotenuse and from the Pythagorean theorem we know that c*c = a*a + b*b. In this case, we are dealing with a segment, usually denoted as s, so a tiny piece of that segment should be called ds. So, the theorem is just ds^2 = dx^2 + dy^2. And the arc length is just the integral of all the tiny segment lengths (ds's).

For those inclined to calculus, this is the integral we would like to compute:

s  =   a sx2
u
x1
V
(1 + (dy/dx)^2) dx
where the limits of integration are calculated as follows:
x1(lat)  =  
a * cos (lat1)
 
V(1 - (2f-f^2)*sin(lat1)^2)
x2(lat2)  =  
a * cos (lat2)
 
V(1 - (2f-f^2)*sin(lat2)^2)

The following code snippet evaluates the above integral numerically. The first step is to calculate the limits of integration, then an appropriate x-increment is calculated, then for each increment the code adds the contribution of Ds = a * Dx * sqrt(1 + (Dy/Dx)^2).

double EllipseArcLength(double lat1, double lat2,
                        double a = GEO::ERAD,
                        double f = GEO::FLATTENING)
{
   double result = 0.0;


   // how many steps to use

   INT32 steps = 100 + 100 * (INT32)(0.50 + (lat2>lat1) ?
                             (lat2-lat1) : (lat1-lat2));
   steps = (steps > 4000.0) ? 4000.0 : steps;

   double snLat1 = sin(GEO::DE2RA*lat1);
   double snLat2 = sin(GEO::DE2RA*lat2);
   double twoF   = 2 * f - f * f;

   // limits of integration
   double x1 = a * cos(GEO::DE2RA * lat1) /
                   sqrt(1 - twoF*snLat1*snLat1);
   double x2 = a * cos(GEO::DE2RA * lat2) /
                   sqrt(1 - twoF*snLat2*snLat2);

   double dx = (x2 - x1) / (double)(steps - 1);
   double x, y1, y2, dy, dydx;
   double adx = (dx < 0.0) ? -dx : dx;    // absolute value of dx

   double a2 = a * a;
   double oneF = 1 - f;

   // now loop through each step adding up all the little
   // hypotenuses
   for (INT32 i = 0; i < (steps - 1); i++){
      x = x1 + dx * i;
      dydx = ((a * oneF * sqrt((1.0 - ((x+dx)*(x+dx))/a2))) -
         (a * oneF * sqrt((1.0 - (x*x)/a2)))) / dx;
      result += adx * sqrt(1.0 + dydx*dydx);
   }

   return result;
}

Normally, when dealing with an ellipse, you start with the lengths of the semi-major and semi-minor axes (a and b). From a and b, you can easily calculate e, and e can easily be transformed into f. The following transformations can be used to go between e, f, a, and b. These transformation would be useful if you wanted to calculate elliptical arc length in terms of e (the eccentricity of an ellipse), rather than f.

f = 1 - (b/a) ; b = a(1 - f) ; e^2 = 2f - f^2 ; f = 1 - sqrt(1-e^2) ; e^2 = sqrt(1 - (b^2/a^2))

Problem 2B. Calculate azimuth and path length (the geodetic forward and inverse problems)

The geodetic forward and inverse problems are the same as Problem 1B and 1A, respectively, only for the ellipsoid instead of a sphere. The basic approach to solving the problems for the azimuth and the distance is similar to Problems 1A and 1B in that triangles on a curved surface are analyzed. The difference is that instead of looking at one single triangle, the path is broken into small (infinitesmal) parts using differential geometry, then by applying theorems from calculus you end up with elliptic integrals that can be solved for the sides or angles of the small curved triangles. Elliptic integrals are usually solved numerically by expanding the differential triangle elements using series approximations. Luckily for us programmers, these problems have been solved and coded. Usually, all we have to do is know how and when to apply the different formulations. The "exact" ellipsoid calculations are included in the sample project.

Code for Approximate Ellipsoid Distance Calculation

double ApproxDistance(double lat1, double lon1, double lat2,
                      double lon2)
{
   lat1 = GEO::DE2RA * lat1;
   lon1 = -GEO::DE2RA * lon1;
   lat2 = GEO::DE2RA * lat2;
   lon2 = -GEO::DE2RA * lon2;

   double F = (lat1 + lat2) / 2.0;
   double G = (lat1 - lat2) / 2.0;
   double L = (lon1 - lon2) / 2.0;

   double sing = sin(G);
   double cosl = cos(L);
   double cosf = cos(F);
   double sinl = sin(L);
   double sinf = sin(F);
   double cosg = cos(G);

   double S = sing*sing*cosl*cosl + cosf*cosf*sinl*sinl;
   double C = cosg*cosg*cosl*cosl + sinf*sinf*sinl*sinl;
   double W = atan2(sqrt(S),sqrt(C));
   double R = 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 * GEO::ERAD;
   return (D * (1 + GEO::FLATTENING * H1 * sinf*sinf*cosg*cosg -
   GEO::FLATTENING*H2*cosf*cosf*sing*sing));
}

Although the preceding code is called approximate, it is actually much more accurate than the great circle calculation.

Sample Project

The approximate method above can be found in a book by Jean Meeus called Astronomical Algorithms, a terrific book for programmers (he used a formula developed by Andoyer in a publication from 1950 called Annuaire du Bureau des Longitudes). I hope this article and the sample project are helpful for developers who wish to make accurate geographic calculations. Also, I want to thank CodeGuru for posting this article and Ted Yezek for taking his valuable time to discuss and explain so many geographic calculations and GIS algorithms. Please keep in mind that the code in this article and in the sample project were developed for the purpose of explaining the calculations; they have not been through the type of rigorous testing and validation process that professional software is subjected to before release. The "intersections" problem was the result of a reader's question, so if you have other geographic or geometric problems you would like included in the article, send me an e-mail.

I used the following references to write this article:

Astronomical Algorithms by Jean Meeus

Explanatory Supplement to the Astronomical Almanac edited by P. Kenneth Seidelmann

Satellie Communications by Dennis Roddy

Enjoy.



About the Author

Andy McGovern

Andy McGovern - Software Developer/Engineer. Special interests: astronomy, image processing, orbital mechanics, mathematics. Currently work at the Johns Hopkins University Applied Physics Laboratory on the science operations team for the CRISM instrument on the Mars Reconnaissance Orbiter.

Downloads

Comments

  • asin range

    Posted by Anrii on 02/07/2014 06:50am

    Hello Andy. Thanks for your article, it's very useful event after 10 years :) However I've mentioned that using 'asin' to calculate longitude in 'GCIntersectSegment' returns values only in [-pi/2, pi/2] range. While the actual range may be [-pi, pi]. I've replaced it with 'atan2' function like this: atan2((v3a[2] / cos(lat3A)), (v3a[0] / cos(lat3A))); Thanks for your great article once again.

    Reply
  • Problem of intersercting two geodesic paths to find new geographic coordinates

    Posted by Jorge Taramona on 02/02/2014 11:21am

    Congratulations Andy for your webpage, it is absolutely amazing. I would like to know if you have some explanation about math and implementations related to find a geographical coordinate by intersecting two geodesic paths in a similar fashion as intersecting two great circles. Thanks a lot in advance.

    Reply
  • thx

    Posted by seby on 03/25/2013 06:57am

    Thanks a lot for this code! however i think there's a tiny mistake in GCIntersectSegment() when you're using asin from math.h in C++, it returns a result in radians instead of degrees so a conversion is needed: double d1a3a = ApproxDistance(lat1A, lon1A, lat3A*GEO::RA2DE, lon3A*GEO::RA2DE); double d1b3a = ApproxDistance(lat1B, lon1B, lat3A*GEO::RA2DE, lon3A*GEO::RA2DE); double d2a3a = ApproxDistance(lat2A, lon2A, lat3A*GEO::RA2DE, lon3A*GEO::RA2DE); double d2b3a = ApproxDistance(lat2B, lon2B, lat3A*GEO::RA2DE, lon3A*GEO::RA2DE); double d1a3b = ApproxDistance(lat1A, lon1A, lat3B*GEO::RA2DE, lon3B*GEO::RA2DE); double d1b3b = ApproxDistance(lat1B, lon1B, lat3B*GEO::RA2DE, lon3B*GEO::RA2DE); double d2a3b = ApproxDistance(lat2A, lon2A, lat3B*GEO::RA2DE, lon3B*GEO::RA2DE); double d2b3b = ApproxDistance(lat2B, lon2B, lat3B*GEO::RA2DE, lon3B*GEO::RA2DE); Also to test if the intersection is within the two path, you also need to compare if the sum of the two distances between the path ends and the the intersection point is equal to the path distances (i added a margin but a more correct way might be to calculate the cross track error ): if ((d1a3a

    • Seby can you comment?

      Posted by H. Smith on 03/26/2013 02:59pm

      Start1 End1 1-12-38 N 0-13-8 N 173-0-8 E 170-39-9 W Start2 End2 4-13-29 S 4-19-41 N 176-14-32 W 178-26-10 W The result is way off. Can you help?

      Reply
    • Question for seby

      Posted by H. Smith on 03/26/2013 02:41pm

      Seby, I made the adjustments that you gave above but I am still having an issue with the intersection routine. I am assuming the format is HH-MM-SS, is that correct? For example if I use the following coordinates as a test: Starting Ending Lat: 1-12-38 N 0-13-8 N Lon: 173-0-8 E 170-39-9 W Lat: 4-13-29 S 4-19-41 N Long:176-14-32 W 178-26-10 W The output is: 0.64063 -2.51012 Can you help?

      Reply
    Reply
  • Calculate 2nd GeoPoint

    Posted by Nixit on 12/12/2012 01:24am

    Hello Andy McGovern, Thanks for wonderful article, it helped me a lot with my project. I got one query regarding calculating 2nd point, I got distance, 1st Geopoint and either latitude or longitude of 2nd geopoint, so how can I get either longitude or latitude of the 2nd geopoint? any hint or advice can also work. and thank you for sharing this article

    Reply
  • Thanks

    Posted by Anonymous on 06/05/2012 09:34pm

    Thanks, very useful!

    Reply
Leave a Comment
  • Your email address will not be published. All fields are required.

Top White Papers and Webcasts

  • The explosion in mobile devices and applications has generated a great deal of interest in APIs. Today's businesses are under increased pressure to make it easy to build apps, supply tools to help developers work more quickly, and deploy operational analytics so they can track users, developers, application performance, and more. Apigee Edge provides comprehensive API delivery tools and both operational and business-level analytics in an integrated platform. It is available as on-premise software or through …

  • The first phase of API management was about realizing the business value of APIs. This next wave of API management enables the hyper-connected enterprise to drive and scale their businesses as API models become more complex and sophisticated. Today, real world product launches begin with an API program and strategy in mind. This API-first approach to development will only continue to increase, driven by an increasingly interconnected web of devices, organizations, and people. To support this rapid growth, …

Most Popular Programming Stories

More for Developers

Latest Developer Headlines

RSS Feeds