8

I have been using the Moveable-Type website to aid me in some Geocoordinate calcuations and it's been very useful, however, I have a bug in my calculation of the mid-point between two coordinates. My result is close to the expected, but not close enough:

posA = {47.64570362, -122.14073746}
posB = {47.64316917, -122.14032175}

expected result (taken from the movable type calculator) = 47°38′40″N, 122°08′26″W = {47.644444, -122.140556} my result: {49.6054801645915, -122.14052959995759}

Here is my code:

private Geocoordinate MidPoint(Geocoordinate posA, Geocoordinate posB)
{
   Geocoordinate midPoint = new Geocoordinate();

   double dLon = DegreesToRadians(posB.Longitude - posA.Longitude);
   double Bx = Math.Cos(DegreesToRadians(posB.Latitude)) * Math.Cos(dLon);
   double By = Math.Cos(DegreesToRadians(posB.Latitude)) * Math.Sin(dLon);

   midPoint.Latitude = RadiansToDegrees(Math.Atan2(Math.Sin(DegreesToRadians(posA.Latitude)) + Math.Sin(DegreesToRadians(posB.Latitude)), 
                Math.Sqrt((Math.Cos(DegreesToRadians(posA.Latitude)) + Bx) * (Math.Cos(DegreesToRadians(posA.Latitude))) + Bx) + By * By));

   midPoint.Longitude = posA.Longitude + RadiansToDegrees(Math.Atan2(By, Math.Cos(DegreesToRadians(posA.Latitude)) + Bx));

   return midPoint;
}

I've got a couple of private methods to do the conversion between Degrees and Radians and back. E.g.

private double DegreeToRadian(double angle)
{
   return Math.PI * angle / 180.0;
}

I can't work out why my results are off by a couple of degrees on the Lat value. Any ideas?

Thanks

Brad
  • 15,361
  • 6
  • 36
  • 57
Stevieboy84
  • 155
  • 1
  • 9
  • Off my **whole** degrees (you expect `34.8954` but get `29.8954`) or **approximate** (you expect `34.8954` but get `29.7836`)? – Brad Nov 12 '10 at 13:20
  • i'm not quite sure if you're considering the earth's surface as a simple sphere or an ellipsoidal sphere. the last time i worked w/ Lambert's law [using a Java implementation] and it seemed to work accurately. – anirvan Nov 12 '10 at 13:26
  • Thanks. It is out by 2 degrees on the lat, and the lon is pretty much spot on. I'm using the formula for Mid-point calculation on this website, translated into C# http://www.movable-type.co.uk/scripts/latlong.html – Stevieboy84 Nov 12 '10 at 14:36
  • I hate trig... Thank you for this! – Rikon Mar 22 '13 at 20:24

2 Answers2

10

You placed some parentheses wrong. I marked the place in the code.

private Geocoordinate MidPoint(Geocoordinate posA, Geocoordinate posB)
{
   Geocoordinate midPoint = new Geocoordinate();

   double dLon = DegreesToRadians(posB.Longitude - posA.Longitude);
   double Bx = Math.Cos(DegreesToRadians(posB.Latitude)) * Math.Cos(dLon);
   double By = Math.Cos(DegreesToRadians(posB.Latitude)) * Math.Sin(dLon);

   midPoint.Latitude = RadiansToDegrees(Math.Atan2(
                Math.Sin(DegreesToRadians(posA.Latitude)) + Math.Sin(DegreesToRadians(posB.Latitude)),
                Math.Sqrt(
                    (Math.Cos(DegreesToRadians(posA.Latitude)) + Bx) *
                    (Math.Cos(DegreesToRadians(posA.Latitude)) + Bx) + By * By))); 
                 // (Math.Cos(DegreesToRadians(posA.Latitude))) + Bx) + By * By)); // Your Code

   midPoint.Longitude = posA.Longitude + RadiansToDegrees(Math.Atan2(By, Math.Cos(DegreesToRadians(posA.Latitude)) + Bx));

   return midPoint;
}
PetPaulsen
  • 3,442
  • 2
  • 22
  • 33
0

Some methods (eg WGS84 conventions) include oblateness of the earth. (at least that's what it says here)

Sanjay Manohar
  • 6,920
  • 3
  • 35
  • 58