Python lat/long midpoint calculation gives wrong result when longitude > 90

Question

I have a problem with a short function to calculate the midpoint of a line when given the latitude and longitude of the points at each end. To put it simply, it works correctly when the longitude is greater than -90 degrees or less than 90 degrees. For the other half of the planet, it provides a somewhat random result.

The code is a python conversion of javascript provided at http://www.movable-type.co.uk/scripts/latlong.html, and appears to conform to the corrected versions here and here. When comparing with the two stackoverflow versions, I'll admit I don't code in C# or Java, but I can't spot where my error is.

Code is as follows:

#!/usr/bin/python

import math

def midpoint(p1, p2):
   lat1, lat2 = math.radians(p1[0]), math.radians(p2[0])
   lon1, lon2 = math.radians(p1[1]), math.radians(p2[1])
   dlon = lon2 - lon1
   dx = math.cos(lat2) * math.cos(dlon)
   dy = math.cos(lat2) * math.sin(dlon)
   lat3 = math.atan2(math.sin(lat1) + math.sin(lat2), math.sqrt((math.cos(lat1) + dx) * (math.cos(lat1) + dx) + dy * dy))
   lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx)
   return(math.degrees(lat3), math.degrees(lon3))

p1 = (6.4, 45)
p2 = (7.3, 43.5)
print "Correct:", midpoint(p1, p2)

p1 = (95.5,41.4)
p2 = (96.3,41.8)
print "Wrong:", midpoint(p1, p2)

Any suggestions?

Take comfort: http://www.ig.utexas.edu/outreach/googleearth/latlong.html is also wrong. — S.Lott, May 05 '11 at 11:09
@John Machin: The second example provides an answer that does not appear to be between the two points. — S.Lott, May 05 '11 at 21:34
@S.Lott: If you mean the OP's second example p1=(95.5,41.4) etc: 95.5 degrees of LATITUDE is invalid -- it would be further north than the North Pole. Read my answer -- the OP got lat and lon reversed and has confessed. If you mean something else, please explain. — John Machin, May 05 '11 at 22:06
@S.Lott: Yes both should raise an exception if the args are invalid. They could also ensure that their output is in the valid range -- see the update to my answer. However this is nothing to do with the "same error for short distances" that you haven't explained yet. — John Machin, May 05 '11 at 23:52
Just to clarify - John is correct; it was a schoolboy error, but as it had provided valid results in all previous usage, I hadn't picked it up. I understand the problem of accuracy over short distances (perfect sphere vs. WGS84), and I do use haversine for distances. However, the error factor on midpoints isn't significant for my purposes. — ichneumonad, May 06 '11 at 00:02
@John Machin: "...that you haven't explained yet". I tried several times to explain. I even provided the example of my mistaken understanding. What more do you want? How many additional examples of my mistaken understanding do you want? — S.Lott, May 06 '11 at 02:57
@S.Lott: Sorry. I didn't understand that what you were saying was an example of your mistaken understanding. — John Machin, May 06 '11 at 04:51
@John Machin: "both should raise an exception if the args are invalid" indicated you did actually understand my mistake. — S.Lott, May 06 '11 at 10:00
@S.Lott: "invalid args" has no correlation with "short distances". In any case, this conversation is all a bit pointless -- shall we delete it? — John Machin, May 06 '11 at 10:47
@John Machin: Correct. And you clearly spotted and corrected my error. While *also* asking what my error was. Confusing. — S.Lott, May 06 '11 at 11:13

John Machin · Accepted Answer · 2011-05-05T23:40:01.677

Replace your arg set up code by:

lat1, lon1 = p1
lat2, lon2 = p2
assert -90 <= lat1 <= 90
assert -90 <= lat2 <= 90
assert -180 <= lon1 <= 180
assert -180 <= lon2 <= 180
lat1, lon1, lat2, lon2 = map(math.radians, (lat1, lon1, lat2, lon2))

and run your code again.

Update A few hopefully-helpful general suggestions about calculations involving latitude/longitude:

Input lat/lon in degrees or radians?
Check input lat/lon for valid range
Check OUTPUT lat/lon for valid range. Longitude has a discontinuity at the international dateline.

The last part of the midpoint routine could be usefully changed to avoid a potential problem with long-distance use:

lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx)
# replacement code follows:
lon3d = math.degrees(lon3)
if lon3d < -180:
    print "oops1", lon3d
    lon3d += 360
elif lon3d > 180:
    print "oops2", lon3d
    lon3d -= 360
return(math.degrees(lat3), lon3d)

For example, finding a midpoint between Auckland, New Zealand (-36.9, 174.8) and Papeete, Tahiti (-17.5, -149.5) produces oops2 194.270430902 on the way to a valid answer (-28.355951246746923, -165.72956909809082)

corrected to assert p1[0] etc. and now I feel like an idiot - yep, the latitude and longitude were the wrong way round. Still provided the right answers for half the planet though! Thanks John. — ichneumonad, May 05 '11 at 10:58
@ichneumonad: refer to my edited version for suggested better style than all that p0[1] stuff — John Machin, May 05 '11 at 11:02
@ichneumonad: I'm finding it a bit hard to believe that swapping lat and lon gives the same answer for "half the planet". For example: `midpoint((1, 2), (3, 4))` produces `(2.0003044085023722, 2.999390393801055)` but `midpoint((2, 1), (4, 3))` produces `(3.0004561487854735, 1.9990851259125342)` — John Machin, May 09 '11 at 00:42

score 1 · Answer 2 · answered May 14 '12 at 15:16

First of all, apologies that I'm leaving another answer. I have a solution to the problem mentioned in that answer involving how to find the midpoint of two points on either side of the dateline. I'd love to simply add a comment to the existing answer, but I don't have the reputation to do so.

The solution was found by looking at the Javascript files powering the tool at http://www.movable-type.co.uk/scripts/latlong.html. I'm using shapely.geometry.Point. If you don't want to install this package, then using tuples instead would work just the same.

def midpoint(pointA, pointB):
    lonA = math.radians(pointA.x)
    lonB = math.radians(pointB.x)
    latA = math.radians(pointA.y)
    latB = math.radians(pointB.y)

    dLon = lonB - lonA

    Bx = math.cos(latB) * math.cos(dLon)
    By = math.cos(latB) * math.sin(dLon)

    latC = math.atan2(math.sin(latA) + math.sin(latB),
                  math.sqrt((math.cos(latA) + Bx) * (math.cos(latA) + Bx) + By * By))
    lonC = lonA + math.atan2(By, math.cos(latA) + Bx)
    lonC = (lonC + 3 * math.pi) % (2 * math.pi) - math.pi

    return Point(math.degrees(lonC), math.degrees(latC))

I hope this is helpful and not regarded as inappropriate seeing as how it is an answer to the question raised in the previous answer.

score 0 · Answer 3 · answered Jun 20 '17 at 10:10

Not straight answer to the >90 question, but it helped in my case so I leave it here just in case it could help others.

I only needed to place the mid point in a map with lat,long in degrees. I used no conversion to radians and it is correctly depicted in the map by using simple euclidean distance. Example function:

def midpoint_euclidean(x1,y1,x2,y2):
    dist_x = abs(x1-x2) / 2.
    dist_y = abs(y1-y2) / 2.
    res_x = x1 - dist_x if x1 > x2 else x2 - dist_x
    res_y = y1 - dist_y if y1 > y2 else y2 - dist_y
    return res_x, res_y

Python lat/long midpoint calculation gives wrong result when longitude > 90

3 Answers3

Linked