# Find longitude for point between two coordinates with pyproj?

I've written a function to get the longitude of a point of known latitude between two known points (might not be working correctly yet), but I feel like there must be a better way of using pyproj to do this.

There isn't a built in method for doing this, but doing this binary search method I'm using seems possibly unnecessary as opposed to just solving for the other half of the coordinate directly.

Any suggestions? Ideally I'd like to generalize this to solve for latitude when the longitude is known, too.

My function thing (maybe not working, will update once I'm sure):

``````import pyproj
import random

def unknownlon(lon1,lat1,lon2,lat2,lat_known):
''' Binary search by distance to find unknown longitude between points

latitude and azimuth are known'''

# Get positive for north movement and negative for south
vert_change = (lat2 - lat1) / abs(lat2 - lat1)
horz_change = (lon2 - lon1) / abs(lon2 - lon1)

# Calculate angle and distance between known points
fwd_azimuth,back_azimuth,dist = g.inv(lon1,lat1,lon2,lon2)

# Random distance between known points
dist_result = random.uniform(0.0, dist)

# Loop while difference between desired lat and result is high
lat_diff = 1.0
while lat_diff > 0.00000001:
lon_result,lat_result,back_azimuth = g.fwd(lon1,lat1,
fwd_azimuth,dist_result)

# Calc difference between desired lat and result
lat_diff = abs(lat_known - lat_result)

# Determine whether to use a distance higher or lower than last
if (lat_result * vert_change) > (lat_known * vert_change):
dist_result = random.uniform(0.0, dist_result)
else:
dist_result = random.uniform(0.0, dist - dist_result)
lon1 = lon_result
lat1 = lat_result

return lon_result, dist_result, lat_diff

if __name__ == '__main__':

lat_known = 65.0

lat1 = 68.0
lon1 = 3.0
lat2 = 63.0
lon2 = -12.0

g = pyproj.Geod(ellps='WGS84') # Use WGS84 ellipsoid

# Calculate an unknown longitude
lon_result, dist_result, lat_diff = unknownlon(lon1,lat1,
lon2,lat2,
lat_known)
``````
• How accurate must the calculation be? What are typical distances between the points? Do you have to handle points near the poles? What do you mean by "solve for latitude" when all latitudes in this problem are given as inputs? Commented Nov 8, 2013 at 20:47
• I would like it to be as accurate as possible; though I suppose 1km would be fine. I suppose it matters how near the poles, but I could definitely need this near 80N, maybe higher. The last I was just explaining that I would like to generalize this, so I can solve for either lat/lon given a lon/lat and two known points. Commented Nov 8, 2013 at 20:53
• One km is easy when your points are only a few kilometers apart :-). On the other hand, that level of accuracy requires different methods when your points could be across the globe. I'm trying to learn whether you can settle for the relatively simple formulas of spherical (or Euclidean 3D) geometry or if you need the more complicated formulas for ellipsoids. Commented Nov 8, 2013 at 20:56
• The `pyproj` package can use standard ellipsoid definitions for it's calculation, which is why I was thinking it would be a good option. I'll add the class/function call to my question for clarity. Thanks! Commented Nov 9, 2013 at 0:20

Your first problem has a straightforward solution...

Given lat1, lon1, lat2, lon2, you can solve the inverse geodesic problem to determine azi1, azi2, the azimuths at the end points. The problem now is: given lat1, lon1, azi1, and lat3 determine lon3. This is the so-called hybrid problem and its solution is described in Section 4 of Algorithms for geodesics (addenda). First of all, note that depending on the value of lon3 there may be 0, 1, or 2 solutions. However, in the case where lat3 lies between lat1 and lat2, there is a single solution. The solution involves converting the latitudes to parametric latitudes and solving the resulting problem on the auxiliary sphere to give the spherical arc distance between points 1 and 3. You can then solve the direct problem with this given spherical arc distance to determin lon3. Unfortunately the pyproj API doesn't provide hooks to solving the direct problem in this way. You can get around this by using the underlying C interface in proj.4 or else by using the native python solution for the geodesic problem. (Unfortunately, I see the solution to this specific problem is rather buried in the python code. You'll need to use GeodesicLine.GenPosition with arcmode = True. I'll provide a friendlier API with the next release which will appear within the next month. In the meantime, see below for the patch.)

The second problem you posed (determining the lat3 given lon3) is Problem 4 of Section 10 of Geodesics on an ellipsoid of revolution (errata). This can be solved using Newton's method.

P.S. The patch to the python geodesic routines to support specifying distance in terms of arc length is given by this patch.

• Excellent answer. Until I have time to look at these options, would you say that my method is adequate (i.e. not flawed in some way)? Thank you. Commented Nov 11, 2013 at 12:28
• A few comments: (1) It'll only work for lat_known between lat1 and lat2. (2) It's a little quirky to pick the new guess randomly; each time you call it you'll get a slightly different answer. (3) Your else clause is probably wrong: shouldn't you update dist? in addition, you accumulate round-off errors by resetting lat1 and lon1 each time.
– cffk
Commented Nov 11, 2013 at 15:58
• P.S. (4) You'll get a division by zero error if lat1 == lat2 (or if lon1 == lon2).
– cffk
Commented Nov 12, 2013 at 1:41

Here's the complete solution for the latitude intersection problem. It assumes that the python version of GeographicLib is installed and that you've applied ArcDirect patch. Running this script (which solves the particular problem posed in the question), I get "Intersect at longitude -7.16007475, distance 561449.408". The code tries to do something reasonable in all cases. Let me know if that's incorrect.

``````import math

def LatitudeIntersect(lat1,azi1,lat2,f):
"""Return arc distance on auxiliary sphere for a geodesic starting at
lat1 with azimuth to its first intersection with latitude lat2.
Flattening of ellipsoid is f and angles are measured in degrees.
Require lat1,lat2 in [-90,90], azi1 in [-180,180], f < 1.  Return nan
if no intersection.  Return 0 for equatorial case (lat1 = lat2 = 0,
azi1 = +/-90).  Else if lat1 == lat2, return 0 if geodesic is
initially heading towards or parallel to the equator; otherwise return
first positive intersection.
"""
degree = math.pi/180; tiny = math.pow(2.0, -511); nan = float("nan")

phi1 = lat1 * degree
sphi1 = math.sin(phi1); cphi1 = 0 if abs(phi1) == 90 else math.cos(phi1)
sbet1, cbet1 = (1-f)*sphi1, cphi1
r = math.hypot(sbet1, cbet1); sbet1, cbet1 = sbet1/r, cbet1/r

phi2 = lat2 * degree
sphi2 = math.sin(phi2); cphi2 = 0 if abs(phi2) == 90 else math.cos(phi2)
sbet2, cbet2 = (1-f)*sphi2, cphi2
r = math.hypot(sbet2, cbet2); sbet2, cbet2 = sbet2/r, cbet2/r

alp1 = azi1 * degree
salp1 = 0 if abs(azi1) == 180 else math.sin(alp1)
calp1 = 0 if abs(azi1) ==  90 else math.cos(alp1)

if sbet1 == 0 and calp1 == 0: # equatorial case
return 0 if lat2 == 0 else nan

salp0, calp0 = salp1*cbet1, math.hypot(calp1, salp1*sbet1)
ssig1, csig1 = sbet1, calp1*cbet1
r = math.hypot(ssig1, csig1); ssig1, csig1 = ssig1/r, csig1/r

if abs(lat2) == abs(lat1):
salp2, calp2 = salp1, abs(calp1)
else:
salp2 = salp0/max(tiny,cbet2)
calp2 = calp0*calp0 - sbet2*sbet2
if calp2 < 0:
return nan
calp2 = math.sqrt(calp2)/cbet2

if lat2 < lat1:
calp2 = -calp2
elif lat2 == lat1:
if lat2 > 0:
calp2 = -calp2
elif lat2 == 0:
calp2 = -calp1

ssig2, csig2 = sbet2, calp2*cbet2
r = math.hypot(ssig2, csig2); ssig2, csig2 = ssig2/r, csig2/r

ssig12 = ssig2*csig1 - csig2*ssig1; csig12 = csig2*csig1 + ssig2*ssig1
a12 = math.atan2(ssig12, csig12)/degree

return a12 if a12 >= 0 else a12+360

from geographiclib.geodesic import Geodesic
lat1,lon1 = 68,3
lat2,lon2 = 63,-12
lat3 = 65
g12=Geodesic.WGS84.Inverse(lat1,lon1,lat2,lon2)
a13=LatitudeIntersect(68,g12['azi1'],65,Geodesic.WGS84._f)
g13=Geodesic.WGS84.ArcDirect(lat1,lon1,g12['azi1'],a13)
print 'Intersect at longitude %.8f, distance %.3f' % (g13['lon2'], g13['s12'])
``````
• Stellar. That longitude certainly looks correct. I'd still like to make the corrections to my script using `pyproj`. I will add an update with any comparisons to be made. Thank you for the help. Commented Nov 18, 2013 at 10:39
• Version 1.34 of the python package for geodesics is now available at pypi.python.org/pypi/geographiclib/1.34 This includes the "ArcDirect patch".
– cffk
Commented Dec 11, 2013 at 19:09