Find longitude for point between two coordinates with pyproj?

Question

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)

Answer 1

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.

Answer 2

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'])