Skip to content
This repository has been archived by the owner on Jan 19, 2021. It is now read-only.

Commit

Permalink
added a bunch of geometry functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@Guymer
Guymer committed Feb 25, 2017
1 parent dc43c28 commit 32a8720
Show file tree
Hide file tree
Showing 6 changed files with 249 additions and 0 deletions.
5 changes: 5 additions & 0 deletions __init__.py
View file
Expand Up @@ -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
Expand Down
25 changes: 25 additions & 0 deletions calc_angle_between_two_locs.py
View file
@@ -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)
77 changes: 77 additions & 0 deletions calc_dist_between_two_locs.py
View file
@@ -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)
77 changes: 77 additions & 0 deletions calc_loc_from_loc_and_bearing_and_dist.py
View file
@@ -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)
26 changes: 26 additions & 0 deletions find_middle_of_great_circle.py
View file
@@ -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)
39 changes: 39 additions & 0 deletions find_point_on_great_circle.py
View file
@@ -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)

0 comments on commit 32a8720

Please sign in to comment.