From 32a872073977f63e79d9a40a6f7ac306eead8b23 Mon Sep 17 00:00:00 2001 From: Thomas Guymer Date: Sat, 25 Feb 2017 08:58:03 +0000 Subject: [PATCH] added a bunch of geometry functions --- __init__.py | 5 ++ calc_angle_between_two_locs.py | 25 ++++++++ calc_dist_between_two_locs.py | 77 +++++++++++++++++++++++ calc_loc_from_loc_and_bearing_and_dist.py | 77 +++++++++++++++++++++++ find_middle_of_great_circle.py | 26 ++++++++ find_point_on_great_circle.py | 39 ++++++++++++ 6 files changed, 249 insertions(+) create mode 100644 calc_angle_between_two_locs.py create mode 100644 calc_dist_between_two_locs.py create mode 100644 calc_loc_from_loc_and_bearing_and_dist.py create mode 100644 find_middle_of_great_circle.py create mode 100644 find_point_on_great_circle.py diff --git a/__init__.py b/__init__.py index d65b3ea..2ed4ec7 100644 --- a/__init__.py +++ b/__init__.py @@ -6,8 +6,13 @@ """ # Load sub-functions ... +from .calc_angle_between_two_locs import calc_angle_between_two_locs +from .calc_dist_between_two_locs import calc_dist_between_two_locs +from .calc_loc_from_loc_and_bearing_and_dist import calc_loc_from_loc_and_bearing_and_dist from .convert_bytes_to_pretty_bytes import convert_bytes_to_pretty_bytes from .does_file_have_RTP_hints import does_file_have_RTP_hints +from .find_middle_of_great_circle import find_middle_of_great_circle +from .find_point_on_great_circle import find_point_on_great_circle from .find_program_version import find_program_version from .generate_password import generate_password from .generate_random_stub import generate_random_stub diff --git a/calc_angle_between_two_locs.py b/calc_angle_between_two_locs.py new file mode 100644 index 0000000..7b0adab --- /dev/null +++ b/calc_angle_between_two_locs.py @@ -0,0 +1,25 @@ +# -*- coding: utf-8 -*- + +def calc_angle_between_two_locs(lon1_deg, lat1_deg, lon2_deg, lat2_deg): + # NOTE: math.sqrt() has been replaced with math.hypot() where possible. + # NOTE: math.atan() has been replaced with math.atan2() where possible. + + # Import modules ... + import math + + # Convert to radians ... + lon1_rad = math.radians(lon1_deg) # [rad] + lat1_rad = math.radians(lat1_deg) # [rad] + lon2_rad = math.radians(lon2_deg) # [rad] + lat2_rad = math.radians(lat2_deg) # [rad] + + # Calculate angle in radians ... + distance_rad = 2.0 * math.asin( + math.hypot( + math.sin((lat1_rad - lat2_rad) / 2.0), + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin((lon1_rad - lon2_rad) / 2.0) + ) + ) # [rad] + + # Return angle ... + return math.degrees(distance_rad) diff --git a/calc_dist_between_two_locs.py b/calc_dist_between_two_locs.py new file mode 100644 index 0000000..8ce11c0 --- /dev/null +++ b/calc_dist_between_two_locs.py @@ -0,0 +1,77 @@ +# -*- coding: utf-8 -*- + +def calc_dist_between_two_locs(lon1_deg, lat1_deg, lon2_deg, lat2_deg): + # NOTE: https://en.wikipedia.org/wiki/Vincenty%27s_formulae + # NOTE: math.sqrt() has been replaced with math.hypot() where possible. + # NOTE: math.atan() has been replaced with math.atan2() where possible. + # NOTE: "lambda" is a reserved word in Python so I use "lam" as my variable name. + + # Import modules ... + import math + + # Convert to radians ... + lon1 = math.radians(lon1_deg) # [rad] + lat1 = math.radians(lat1_deg) # [rad] + lon2 = math.radians(lon2_deg) # [rad] + lat2 = math.radians(lat2_deg) # [rad] + + # Set constants ... + a = 6378137.0 # [m] + f = 1.0 / 298.257223563 + b = (1.0 - f) * a # [m] + l = lon2 - lon1 # [rad] + u1 = math.atan((1.0 - f) * math.tan(lat1)) # [rad] + u2 = math.atan((1.0 - f) * math.tan(lat2)) # [rad] + + # Set initial value of lambda and initialize counter ... + lam = l + i = 0 + + # Start infinite loop ... + while True: + # Stop looping if the function has been called too many times ... + if i >= 1000: + break + + # Calculate new lambda and increment counter ... + sin_sigma = math.hypot( + math.cos(u2) * math.sin(lam), + math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lam) + ) + cos_sigma = math.sin(u1) * math.sin(u2) + math.cos(u1) * math.cos(u2) * math.cos(lam) + sigma = math.atan2( + sin_sigma, + cos_sigma + ) + sin_alpha = math.cos(u1) * math.cos(u2) * math.sin(lam) / sin_sigma + cosSq_alpha = 1.0 - sin_alpha ** 2 + cos_two_sigma_m = cos_sigma - 2.0 * math.sin(u1) * math.sin(u2) / cosSq_alpha + c = f * cosSq_alpha * (4.0 + f * (4.0 - 3.0 * cosSq_alpha)) / 16.0 + 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))) + i += 1 + + # Only check the solution after at least 3 function calls ... + if i >= 3: + if abs(lamNew - lam) / abs(lamNew) <= 1.0e-12: + break + + # Replace old lambda with new lambda ... + lam = lamNew + + # Calculate ellipsoidal distance and forward azimuths ... + uSq = cosSq_alpha * (a ** 2 - b ** 2) / b ** 2 + bigA = 1.0 + uSq * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))) / 16384.0 + bigB = uSq * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))) / 1024.0 + 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)) + s = b * bigA * (sigma - delta_sigma) + alpha1 = math.atan2( + math.cos(u2) * math.sin(lam), + math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lam) + ) + alpha2 = math.atan2( + math.cos(u1) * math.sin(lam), + math.sin(u1) * math.cos(u2) - math.cos(u1) * math.sin(u2) * math.cos(lam) + ) + + # Return distance and forward azimuths ... + return s, math.degrees(alpha1), math.degrees(alpha2) diff --git a/calc_loc_from_loc_and_bearing_and_dist.py b/calc_loc_from_loc_and_bearing_and_dist.py new file mode 100644 index 0000000..1d04560 --- /dev/null +++ b/calc_loc_from_loc_and_bearing_and_dist.py @@ -0,0 +1,77 @@ +# -*- coding: utf-8 -*- + +def calc_loc_from_loc_and_bearing_and_dist(lon1_deg, lat1_deg, alpha1_deg, s_m): + # NOTE: https://en.wikipedia.org/wiki/Vincenty%27s_formulae + # NOTE: math.sqrt() has been replaced with math.hypot() where possible. + # NOTE: math.atan() has been replaced with math.atan2() where possible. + # NOTE: "lambda" is a reserved word in Python so I use "lam" as my variable name. + + # Import modules ... + import math + + # Convert to radians ... + lon1 = math.radians(lon1_deg) # [rad] + lat1 = math.radians(lat1_deg) # [rad] + alpha1 = math.radians(alpha1_deg) # [rad] + + # Set constants ... + a = 6378137.0 # [m] + f = 1.0 / 298.257223563 + b = (1.0 - f) * a # [m] + u1 = math.atan((1.0 - f) * math.tan(lat1)) # [rad] + sigma1 = math.atan2( + math.tan(u1), + math.cos(alpha1) + ) + sin_alpha = math.cos(u1) * math.sin(alpha1) + cosSq_alpha = 1.0 - sin_alpha ** 2 + c = f * cosSq_alpha * (4.0 + f * (4.0 - 3.0 * cosSq_alpha)) / 16.0 + uSq = cosSq_alpha * (a ** 2 - b ** 2) / b ** 2 + bigA = 1.0 + uSq * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))) / 16384.0 + bigB = uSq * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))) / 1024.0 + + # Set initial value of sigma and initialize counter ... + sigma = s_m / (b * bigA) + i = 0 + + # Start infinite loop ... + while True: + # Stop looping if the function has been called too many times ... + if i >= 1000: + break + + # Find new value of sigma and increment counter ... + two_sigma_m = 2 * sigma1 + sigma + delta_sigma = bigB * math.sin(sigma) * (math.cos(two_sigma_m) + 0.25 * bigB * (math.cos(sigma) * (2.0 * math.cos(two_sigma_m) ** 2 - 1.0) - bigB * math.cos(two_sigma_m) * (4.0 * math.sin(sigma) ** 2 - 3.0) * (4.0 * math.cos(two_sigma_m) ** 2 - 3.0) / 6.0)) + sigmaNew = s_m / (b * bigA) + delta_sigma + i += 1 + + # Only check the solution after at least 3 function calls ... + if i >= 3: + if abs(sigmaNew - sigma) / abs(sigmaNew) <= 1.0e-12: + break + + # Replace old sigma with new sigma ... + sigma = sigmaNew + + # Calculate end point and forward azimuth ... + lat2 = math.atan2( + math.sin(u1) * math.cos(sigma) + math.cos(u1) * math.sin(sigma) * math.cos(alpha1), + (1.0 - f) * math.hypot( + sin_alpha, + math.sin(u1) * math.sin(sigma) - math.cos(u1) * math.cos(sigma) * math.cos(alpha1) + ) + ) + lam = math.atan2( + math.sin(sigma) * math.sin(alpha1), + math.cos(u1) * math.cos(sigma) - math.sin(u1) * math.sin(sigma) * math.cos(alpha1) + ) + l = lam - (1.0 - c) * f * sin_alpha * (sigma + c * math.sin(sigma) * (math.cos(two_sigma_m) + c * math.cos(sigma) * (2.0 * math.cos(two_sigma_m) ** 2 - 1.0))) + lon2 = l + lon1 + alpha2 = math.atan2( + sin_alpha, + math.cos(u1) * math.cos(sigma) * math.cos(alpha1) - math.sin(u1) * math.sin(sigma) + ) + + # Return end point and forward azimuth ... + return math.degrees(lon2), math.degrees(lat2), math.degrees(alpha2) diff --git a/find_middle_of_great_circle.py b/find_middle_of_great_circle.py new file mode 100644 index 0000000..13bb1c9 --- /dev/null +++ b/find_middle_of_great_circle.py @@ -0,0 +1,26 @@ +# -*- coding: utf-8 -*- + +def find_middle_of_great_circle(lon1_deg = 0.0, lat1_deg = 0.0, lon2_deg = 0.0, lat2_deg = 0.0): + # NOTE: math.sqrt() has been replaced with math.hypot() where possible. + # NOTE: math.atan() has been replaced with math.atan2() where possible. + + # Import modules ... + import math + + # Convert to radians ... + lon1_rad = math.radians(lon1_deg) # [rad] + lat1_rad = math.radians(lat1_deg) # [rad] + lon2_rad = math.radians(lon2_deg) # [rad] + lat2_rad = math.radians(lat2_deg) # [rad] + + # Calculate mid-point ... + Bx = math.cos(lat2_rad) * math.cos(lon2_rad - lon1_rad) + By = math.cos(lat2_rad) * math.sin(lon2_rad - lon1_rad) + lat3_rad = math.atan2( + math.sin(lat1_rad) + math.sin(lat2_rad), + math.sqrt((math.cos(lat1_rad) + Bx) * (math.cos(lat1_rad) + Bx) + By * By) + ) # [rad] + lon3_rad = lon1_rad + math.atan2(By, math.cos(lat1_rad) + Bx) # [rad] + + # Return middle point ... + return math.degrees(lon3_rad), math.degrees(lat3_rad) diff --git a/find_point_on_great_circle.py b/find_point_on_great_circle.py new file mode 100644 index 0000000..ef770ac --- /dev/null +++ b/find_point_on_great_circle.py @@ -0,0 +1,39 @@ +# -*- coding: utf-8 -*- + +def find_point_on_great_circle(frac = 0.0, lon1_deg = 0.0, lat1_deg = 0.0, lon2_deg = 0.0, lat2_deg = 0.0): + # NOTE: math.sqrt() has been replaced with math.hypot() where possible. + # NOTE: math.atan() has been replaced with math.atan2() where possible. + + # Import modules ... + import math + + # Load sub-functions ... + from .calc_angle_between_two_locs import calc_angle_between_two_locs + + # Convert to radians ... + lon1_rad = math.radians(lon1_deg) # [rad] + lat1_rad = math.radians(lat1_deg) # [rad] + lon2_rad = math.radians(lon2_deg) # [rad] + lat2_rad = math.radians(lat2_deg) # [rad] + + # Calculate mid-point ... + rad = math.radians(calc_angle_between_two_locs(lon1_deg, lat1_deg, lon2_deg, lat2_deg)) # [rad] + a = math.sin((1.0 - frac) * rad) / math.sin(rad) + b = math.sin(frac * rad) / math.sin(rad) + x = a * math.cos(lat1_rad) * math.cos(lon1_rad) + b * math.cos(lat2_rad) * math.cos(lon2_rad) + y = a * math.cos(lat1_rad) * math.sin(lon1_rad) + b * math.cos(lat2_rad) * math.sin(lon2_rad) + z = a * math.sin(lat1_rad) + b * math.sin(lat2_rad) + lat3_rad = math.atan2( + z, + math.hypot( + x, + y + ) + ) # [rad] + lon3_rad = math.atan2( + y, + x + ) # [rad] + + # Return point point ... + return math.degrees(lon3_rad), math.degrees(lat3_rad)