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

1#!/usr/bin/env python3 

2 

3# Define function ... 

4def find_middle_of_great_circle( 

5 lon1_deg, 

6 lat1_deg, 

7 lon2_deg, 

8 lat2_deg, 

9 /, 

10): 

11 """Calculate the middle of the great circle that connects two coordinates. 

12 

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

14 Earth and calculates the middle of the great circle that connects them, 

15 correctly handling crossing over the anti-meridian (should it occur). 

16 

17 Parameters 

18 ---------- 

19 lon1 : float 

20 the longitude of the first coordinate (in degrees) 

21 lat1 : float 

22 the latitude of the first coordinate (in degrees) 

23 lon2 : float 

24 the longitude of the second coordinate (in degrees) 

25 lat2 : float 

26 the latitude of the second coordinate (in degrees) 

27 

28 Returns 

29 ------- 

30 lon3 : float 

31 the longitude of the middle of the great circle (in degrees) 

32 lat3 : float 

33 the latitude of the middle of the great circle (in degrees) 

34 

35 Notes 

36 ----- 

37 Copyright 2017 Thomas Guymer [1]_ 

38 

39 References 

40 ---------- 

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

42 """ 

43 

44 # Import standard modules ... 

45 import math 

46 

47 # ************************************************************************** 

48 

49 # Check arguments ... 

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

51 return lon1_deg, lat1_deg 

52 

53 # Convert to radians ... 

54 lon1_rad = math.radians(lon1_deg) # [rad] 

55 lat1_rad = math.radians(lat1_deg) # [rad] 

56 lon2_rad = math.radians(lon2_deg) # [rad] 

57 lat2_rad = math.radians(lat2_deg) # [rad] 

58 

59 # Calculate mid-point ... 

60 Bx = math.cos(lat2_rad) * math.cos(lon2_rad - lon1_rad) 

61 By = math.cos(lat2_rad) * math.sin(lon2_rad - lon1_rad) 

62 lat3_rad = math.atan2( 

63 math.sin(lat1_rad) + math.sin(lat2_rad), 

64 math.sqrt((math.cos(lat1_rad) + Bx) * (math.cos(lat1_rad) + Bx) + By * By) 

65 ) # [rad] 

66 lon3_rad = lon1_rad + math.atan2(By, math.cos(lat1_rad) + Bx) # [rad] 

67 

68 # Return mid-point ... 

69 return math.degrees(lon3_rad), math.degrees(lat3_rad)