1

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? – whuber Nov 8 '13 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. – ryanjdillon Nov 8 '13 at 20:53
  • 1
    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. – whuber Nov 8 '13 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! – ryanjdillon Nov 9 '13 at 0:20
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'])
  • 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. – ryanjdillon Nov 18 '13 at 10:39
  • 1
    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 Dec 11 '13 at 19:09
3

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. – ryanjdillon Nov 11 '13 at 12:28
  • 1
    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 Nov 11 '13 at 15:58
  • P.S. (4) You'll get a division by zero error if lat1 == lat2 (or if lon1 == lon2). – cffk Nov 12 '13 at 1:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.