Geek Answers Handbook

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

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

The code is a python-based javascript transformation presented at http://www.movable-type.co.uk/scripts/latlong.html and seems to correspond to the patched versions here and here . Compared to the two versions of stackoverflow, I agree that I am not C # or Java code, but I cannot determine where my error is.

The 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?

+4
source share
3 answers

Replace your arg installation code with:

 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 encouraging and helpful guidelines for latitude / longitude calculations:

  • Lat / lon input in degrees or radian?
  • Check lat / lon input for valid range
  • Check the OUTPUT lat / lon value for a valid range. Longitude has a gap on the international line.

The last part of the midpoint subroutine can be usefully modified to avoid a potential long-distance usage problem:

 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, looking for an intermediate point between Auckland, New Zealand (-36.9, 174.8) and Papeete, Tahiti (-17.5, -149.5) gives oops2 194.270430902 on the way to a valid answer (-28.355951246746923, -165.72956909809082)

+4
source

First of all, apologies that I am leaving another answer. I have a solution to the problem mentioned in this answer on how to find the midpoint of two points on either side of the date. I would just like to add a comment to the existing answer, but I do not have such a reputation.

The solution was found by looking at the Javascript files, including the http://www.movable-type.co.uk/scripts/latlong.html tool. I am using shapely.geometry.Point. If you do not want to install this package, then using tuples instead will work 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 useful and not considered inappropriate, seeing how this is the answer to the question raised in the previous answer.

+1
source

Not a direct answer to questions> 90, but it helped in my case, so I leave it here just in case, if it can help others.

I only needed to place the midpoint on the map with a lat, in degrees in length. I did not use the conversion in radians and is correctly displayed on the map using a simple Euclidean distance. Function example:

 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
0
source

More articles:


All Articles