How to use geographiclib - 10 common examples

"""Solve geodesic problems"""

  GEOGRAPHICLIB_GEODESIC_ORDER = 6
  nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nA3x_ = nA3_
  nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC3x_ = (nC3_ * (nC3_ - 1)) // 2
  nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC4x_ = (nC4_ * (nC4_ + 1)) // 2
  maxit1_ = 20
  maxit2_ = maxit1_ + Math.digits + 10

  tiny_ = math.sqrt(Math.minval)
  tol0_ = Math.epsilon
  tol1_ = 200 * tol0_
  tol2_ = math.sqrt(tol0_)
  tolb_ = tol0_ * tol2_
  xthresh_ = 1000 * tol2_

  CAP_NONE = GeodesicCapability.CAP_NONE
  CAP_C1   = GeodesicCapability.CAP_C1
  CAP_C1p  = GeodesicCapability.CAP_C1p
  CAP_C2   = GeodesicCapability.CAP_C2
  CAP_C3   = GeodesicCapability.CAP_C3
  CAP_C4   = GeodesicCapability.CAP_C4
  CAP_ALL  = GeodesicCapability.CAP_ALL
  CAP_MASK = GeodesicCapability.CAP_MASK

def _C1f(eps, c):
    """Private: return C1."""
    coeff = [
      -1, 6, -16, 32,
      -9, 64, -128, 2048,
      9, -16, 768,
      3, -5, 512,
      -7, 1280,
      -7, 2048,
    ]
    eps2 = Math.sq(eps)
    d = eps
    o = 0
    for l in range(1, Geodesic.nC1_ + 1): # l is index of C1p[l]
      m = (Geodesic.nC1_ - l) // 2        # order of polynomial in eps^2
      c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
      o += m + 2
      d *= eps
  _C1f = staticmethod(_C1f)

def _C2f(eps, c):
    """Private: return C2"""
    coeff = [
      1, 2, 16, 32,
      35, 64, 384, 2048,
      15, 80, 768,
      7, 35, 512,
      63, 1280,
      77, 2048,
    ]
    eps2 = Math.sq(eps)
    d = eps
    o = 0
    for l in range(1, Geodesic.nC2_ + 1): # l is index of C2[l]
      m = (Geodesic.nC2_ - l) // 2        # order of polynomial in eps^2
      c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
      o += m + 2
      d *= eps
  _C2f = staticmethod(_C2f)

def _transit(lon1, lon2):
    """Count crossings of prime meridian for AddPoint."""
    # Return 1 or -1 if crossing prime meridian in east or west direction.
    # Otherwise return zero.
    # Compute lon12 the same way as Geodesic::Inverse.
    lon1 = Math.AngNormalize(lon1)
    lon2 = Math.AngNormalize(lon2)
    lon12, _ = Math.AngDiff(lon1, lon2)
    cross = (1 if lon1 <= 0 and lon2 > 0 and lon12 > 0
             else (-1 if lon2 <= 0 and lon1 > 0 and lon12 < 0 else 0))
    return cross
  _transit = staticmethod(_transit)

def transit(lon1, lon2):
    """Count crossings of prime meridian."""
    # Return 1 or -1 if crossing prime meridian in east or west direction.
    # Otherwise return zero.
    # Compute lon12 the same way as Geodesic::Inverse.
    lon1 = Math.AngNormalize(lon1)
    lon2 = Math.AngNormalize(lon2)
    lon12 = Math.AngDiff(lon1, lon2)
    cross = (1 if lon1 < 0 and lon2 >= 0 and lon12 > 0
             else (-1 if lon2 < 0 and lon1 >= 0 and lon12 < 0 else 0))
    return cross
  transit = staticmethod(transit)

def _transit(lon1, lon2):
    """Count crossings of prime meridian for AddPoint."""
    # Return 1 or -1 if crossing prime meridian in east or west direction.
    # Otherwise return zero.
    # Compute lon12 the same way as Geodesic::Inverse.
    lon1 = Math.AngNormalize(lon1)
    lon2 = Math.AngNormalize(lon2)
    lon12, _ = Math.AngDiff(lon1, lon2)
    cross = (1 if lon1 <= 0 and lon2 > 0 and lon12 > 0
             else (-1 if lon2 <= 0 and lon1 > 0 and lon12 < 0 else 0))
    return cross
  _transit = staticmethod(_transit)

Azimuth =   Calculus['azi2']                 
                        SpeedMeterSecond = DistanceBetweenPoint             #GPS refresh rate is actually 1, change parameter for different rates
                       # Time = 1                                            
                        if RemainToUseMeterTotal == 0:
                            if DistanceBetweenPoint >= meters:
                                decimalSecondToAdd = meters / DistanceBetweenPoint
                                RemainToUseMeter = DistanceBetweenPoint - meters
                                if os.name == 'nt':
                                    t = threading.Thread(target= extract,args=(ffmpeg,i,decimalSecondToAdd,self.videoFile,Directory))
                                    t.start()
                                    
                                else:
                                    os.system(ffmpeg+ ' -ss '+ str(i + decimalSecondToAdd) +
                                              ' -i '+ str(self.videoFile) + ' -frames:v 1 ' + Directory + '_sec_' + str(i) + str(decimalSecondToAdd)[1:4] +'.png')
                                
                                CalculusDirect = Geodesic.WGS84.Direct(latitude1, longitude1, Azimuth,decimalSecondToAdd* SpeedMeterSecond) 
                                X,Y,quotainutile = self.transform_wgs84_to_utm(CalculusDirect['lon2'],CalculusDirect['lat2'] )  
                                Z = ele1 + decimalSecondToAdd*(ele2 - ele1)
                                txtGPSFile.write(str(Directory.split('/')[-1]) + '_sec_'  + str(i) + str(decimalSecondToAdd)[1:4]+'.png,' + ' ' + str(X) + ', ' + str(Y) + ', ' + str(Z) + '\n')
                                while RemainToUseMeter > meters:
                                    decimalSecondToAddMore = meters / SpeedMeterSecond
                                    RemainToUseMeter = RemainToUseMeter - meters
                                    decimalSecondToAdd = decimalSecondToAdd + decimalSecondToAddMore
                                    if os.name == 'nt':
                                        os.popen(ffmpeg + ' -ss '+ str(i + decimalSecondToAdd) +
                                             ' -i '+ str('"'+self.videoFile+'"') + ' -frames:v 1 ' +'"'+ Directory + '_sec_' + str(i) + str(decimalSecondToAdd)[1:4] +'.png'+'"')
                                    else:
                                        os.system(ffmpeg+ ' -ss '+ str(i + decimalSecondToAdd) +
                                                  ' -i '+ str(self.videoFile) + ' -frames:v 1 ' + Directory + '_sec_' + str(i) + str(decimalSecondToAdd)[1:4] +'.png')
                                    
                                    CalculusDirect = Geodesic.WGS84.Direct(latitude1, longitude1, Azimuth,decimalSecondToAdd* SpeedMeterSecond) 
                                    X,Y,quotainutile = self.transform_wgs84_to_utm(CalculusDirect['lon2'],CalculusDirect['lat2'] )

def TestEdge(self, azi, s, reverse, sign):
    """Return the results for a tentative additional edge."""
    if self._num == 0:               # we don't have a starting point!
      return 0, Math.nan, Math.nan
    num = self._num + 1
    perimeter = self._perimetersum.Sum() + s
    if self._polyline:
      return num, perimeter, Math.nan

    tempsum =  self._areasum.Sum()
    crossings = self._crossings
    _, lat, lon, _, _, _, _, _, S12 = self._earth.GenDirect(
      self._lat1, self._lon1, azi, False, s, self._mask)
    tempsum += S12
    crossings += PolygonArea.transit(self._lon1, lon)
    _, s12, _, _, _, _, _, S12 = self._earth.GenInverse(
      lat, lon, self._lat0, self._lon0, self._mask)
    perimeter += s12
    tempsum += S12
    crossings += PolygonArea.transit(lon, self._lon0)

def _transit(lon1, lon2):
    """Count crossings of prime meridian for AddPoint."""
    # Return 1 or -1 if crossing prime meridian in east or west direction.
    # Otherwise return zero.
    # Compute lon12 the same way as Geodesic::Inverse.
    lon1 = Math.AngNormalize(lon1)
    lon2 = Math.AngNormalize(lon2)
    lon12, _ = Math.AngDiff(lon1, lon2)
    cross = (1 if lon1 <= 0 and lon2 > 0 and lon12 > 0
             else (-1 if lon2 <= 0 and lon1 > 0 and lon12 < 0 else 0))
    return cross
  _transit = staticmethod(_transit)

dummy, m12b, m0, dummy, dummy = self.Lengths(
          self._n, math.pi + bet12a, sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
          cbet1, cbet2, False, C1a, C2a)
        x = -1 + m12b / (cbet1 * cbet2 * m0 * math.pi)
        betscale = (sbet12a / x if x < -0.01
                    else -self._f * Math.sq(cbet1) * math.pi)
        lamscale = betscale / cbet1
        y = (lam12 - math.pi) / lamscale

      if y > -Geodesic.tol1_ and x > -1 - Geodesic.xthresh_:
        # strip near cut
        if self._f >= 0:
          salp1 = min(1.0, -x); calp1 = - math.sqrt(1 - Math.sq(salp1))
        else:
          calp1 = max((0.0 if x > -Geodesic.tol1_ else -1.0), x)
          salp1 = math.sqrt(1 - Math.sq(calp1))
      else:
        # Estimate alp1, by solving the astroid problem.
        #
        # Could estimate alpha1 = theta + pi/2, directly, i.e.,
        #   calp1 = y/k; salp1 = -x/(1+k);  for _f >= 0
        #   calp1 = x/(1+k); salp1 = -y/k;  for _f < 0 (need to check)
        #
        # However, it's better to estimate omg12 from astroid and use
        # spherical formula to compute alp1.  This reduces the mean number of
        # Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
        # (min 0 max 5).  The changes in the number of iterations are as
        # follows:
        #
        # change percent
        #    1       5
        #    0      78