Coverage for pyguymer3/geo/calc_dist_between_two_locs.py: 92%

1#!/usr/bin/env python3 

2 

3# Define function ... 

4def calc_dist_between_two_locs( 

5 lon1_deg, 

6 lat1_deg, 

7 lon2_deg, 

8 lat2_deg, 

9 /, 

10 *, 

11 eps = 1.0e-12, 

12 nIter = 100, 

13): 

14 """Calculate the distance between two coordinates. 

15 

16 This function reads in two coordinates (in degrees) on the surface of the 

17 Earth and calculates the Geodesic distance (in metres) between them and the 

18 headings (in degrees) from each coordinate to the other one. 

19 

20 Parameters 

21 ---------- 

22 lon1_deg : float 

23 the longitude of the first coordinate (in degrees) 

24 lat1_deg : float 

25 the latitude of the first coordinate (in degrees) 

26 lon2_deg : float 

27 the longitude of the second coordinate (in degrees) 

28 lat2_deg : float 

29 the latitude of the second coordinate (in degrees) 

30 eps : float, optional 

31 the tolerance of the Vincenty formula iterations 

32 nIter : int, optional 

33 the maximum number of iterations (particularly the Vincenty formula) 

34 

35 Returns 

36 ------- 

37 s_m : float 

38 the distance between the two coordinates (in metres) 

39 alpha1_deg : float 

40 the heading to the second coordinate from the first coordinate (in 

41 degrees) 

42 alpha2_deg : float 

43 the heading to the first coordinate from the second coordinate (in 

44 degrees) 

45 

46 Notes 

47 ----- 

48 This function uses `Vincenty's formulae 

49 <https://en.wikipedia.org/wiki/Vincenty%27s_formulae>`_ ; there is a 

50 `JavaScript implementation 

51 <https://www.movable-type.co.uk/scripts/latlong-vincenty.html>`_ online too. 

52 

53 ``lambda`` is a reserved word in Python so I use ``lam`` as my variable name 

54 instead. 

55 

56 Copyright 2017 Thomas Guymer [1]_ 

57 

58 References 

59 ---------- 

60 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3 

61 """ 

62 

63 # Import standard modules ... 

64 import math 

65 

66 # ************************************************************************** 

67 

68 # Skip if the start- and end-points are the same ... 

69 if lon1_deg == lon2_deg and lat1_deg == lat2_deg: 

70 return 0.0, 0.0, 0.0 

71 

72 # Convert to radians ... 

73 lon1 = math.radians(lon1_deg) # [rad] 

74 lat1 = math.radians(lat1_deg) # [rad] 

75 lon2 = math.radians(lon2_deg) # [rad] 

76 lat2 = math.radians(lat2_deg) # [rad] 

77 

78 # Set constants ... 

79 a = 6378137.0 # [m] 

80 f = 1.0 / 298.257223563 

81 b = (1.0 - f) * a # [m] 

82 l = lon2 - lon1 # [rad] 

83 u1 = math.atan((1.0 - f) * math.tan(lat1)) # [rad] 

84 u2 = math.atan((1.0 - f) * math.tan(lat2)) # [rad] 

85 

86 # Set initial value of lambda and initialize counter ... 

87 lam = l 

88 iIter = 0 # [#] 

89 

90 # Start infinite loop ... 

91 while True: 

92 # Stop looping if the function has been called too many times ... 

93 if iIter >= nIter: 

94 raise Exception(f"failed to converge; loc1 = ({lon1_deg:+.9f}°,{lat1_deg:+.9f}°); loc2 = ({lon2_deg:+.9f}°,{lat2_deg:+.9f}°); eps = {eps:.15e}; nIter = {nIter:,d}") from None 

95 

96 # Calculate new lambda and increment counter ... 

97 sin_sigma = math.hypot( 

98 math.cos(u2) * math.sin(lam), 

99 math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lam) 

100 ) 

101 if sin_sigma == 0.0: 

102 # NOTE: co-incident points 

103 return 0.0, 0.0, 0.0 

104 cos_sigma = math.sin(u1) * math.sin(u2) + math.cos(u1) * math.cos(u2) * math.cos(lam) 

105 sigma = math.atan2( 

106 sin_sigma, 

107 cos_sigma 

108 ) 

109 sin_alpha = math.cos(u1) * math.cos(u2) * math.sin(lam) / sin_sigma 

110 cosSq_alpha = 1.0 - sin_alpha ** 2 

111 cos_two_sigma_m = cos_sigma - 2.0 * math.sin(u1) * math.sin(u2) / cosSq_alpha 

112 if math.isnan(cos_two_sigma_m): 

113 # NOTE: equatorial line 

114 cos_two_sigma_m = 0.0 

115 c = f * cosSq_alpha * (4.0 + f * (4.0 - 3.0 * cosSq_alpha)) / 16.0 

116 lamNew = l + (1.0 - c) * f * sin_alpha * (sigma + c * sin_sigma * (cos_two_sigma_m + c * cos_sigma * (2.0 * cos_two_sigma_m ** 2 - 1.0))) 

117 iIter += 1 # [#] 

118 

119 # Only check the solution after at least 3 function calls ... 

120 if iIter >= 3: 

121 if lamNew == lam: 

122 # Replace old lambda with new lambda ... 

123 lam = lamNew 

124 

125 break 

126 if abs(lamNew - lam) / abs(lamNew) <= eps: 

127 # Replace old lambda with new lambda ... 

128 lam = lamNew 

129 

130 break 

131 

132 # Replace old lambda with new lambda ... 

133 lam = lamNew 

134 

135 # Calculate ellipsoidal distance and forward azimuths ... 

136 uSq = cosSq_alpha * (a ** 2 - b ** 2) / b ** 2 

137 bigA = 1.0 + uSq * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))) / 16384.0 

138 bigB = uSq * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))) / 1024.0 

139 delta_sigma = bigB * sin_sigma * (cos_two_sigma_m + 0.25 * bigB * (cos_sigma * (2.0 * cos_two_sigma_m ** 2 - 1.0) - bigB * cos_two_sigma_m * (4.0 * sin_sigma ** 2 - 3.0) * (4.0 * cos_two_sigma_m ** 2 - 3.0) / 6.0)) 

140 s = b * bigA * (sigma - delta_sigma) # [m] 

141 alpha1 = math.atan2( 

142 math.cos(u2) * math.sin(lam), 

143 math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lam) 

144 ) # [rad] 

145 alpha2 = math.atan2( 

146 math.cos(u1) * math.sin(lam), 

147 math.sin(u1) * math.cos(u2) - math.cos(u1) * math.sin(u2) * math.cos(lam) 

148 ) # [rad] 

149 alpha1 = (alpha1 + 2.0 * math.pi) % (2.0 * math.pi) # NOTE: Normalize to +0° <--> +360° 

150 alpha2 = (alpha2 + 2.0 * math.pi) % (2.0 * math.pi) # NOTE: Normalize to +0° <--> +360° 

151 

152 # Return distance and forward azimuths ... 

153 return s, math.degrees(alpha1), math.degrees(alpha2)