Geodesic Distance

distance_with_flattening.py

#!/usr/bin/env python3
import math
from geopy.geocoders import Nominatim
from geopy.distance import great_circle
from geopy.distance import geodesic
from geopy.distance import vincenty
geolocator = Nominatim( user_agent="lab6" )

# https://stackoverflow.com/a/52371976
def decimal_degrees_to_dms( deg , type='lat' ):
	decimals, number = math.modf(deg)
	d = int(number)
	m = int(decimals * 60)
	s = (deg - d - m / 60) * 3600.00
	compass = {
		'lat': ('North','South'),
		'lon': ('East','West')
	}
	compass_str = compass[type][0 if d >= 0 else 1]
	return f"{abs(d)}:{abs(m)}:{abs(s)} {compass_str}"

# https://github.com/kamtechie/DirectionFinder/blob/master/DirectionFinder.py#L28
def compass_direction( lat1 , lon1 , lat2 , lon2 ):
	radians = math.atan2( abs( lon2 - lon1 ) , abs( lat2 - lat1 ) )
	compass_reading = radians * ( 180 / math.pi )
	return compass_reading

# https://www.rdocumentation.org/packages/geosphere/versions/1.5-10/topics/distGeo
def compute_distance_along_meridian( lat1 , lon1 , elevation_1_in_km , lat2 , lon2 , elevation_2_in_km ):
	earth_radius_in_km = 6371
	WGS84_earth_major_equatorial_radius_of_ellipsoid = 6378137
	WGS84_earth_ellipsoid_flattening = float( 1.0 / 298.257223563 )

	lat1_radians = lat1 * ( math.pi/180.0 )
	lon1_radians = lon1 * ( math.pi/180.0 )
	lat2_radians = lat2 * ( math.pi/180.0 )
	lon2_radians = lon2 * ( math.pi/180.0 )

	sin_of_lat1 = math.sin( lat1_radians )
	sin_of_lat2 = math.sin( lat2_radians )

	integration_steps = 40000

	two_f = ( 2 * WGS84_earth_ellipsoid_flattening ) - ( WGS84_earth_ellipsoid_flattening * WGS84_earth_ellipsoid_flattening )

	# Integration Limits
	x1 = ( WGS84_earth_major_equatorial_radius_of_ellipsoid * math.cos( lat1_radians ) ) / math.sqrt( 1 - two_f * sin_of_lat1 * sin_of_lat1 )
	x2 = ( WGS84_earth_major_equatorial_radius_of_ellipsoid * math.cos( lat2_radians ) ) / math.sqrt( 1 - two_f * sin_of_lat2 * sin_of_lat2 )

	dx = ( x2 - x1 ) / ( integration_steps - 1 )
	adx = abs( dx )

	a_squared = ( WGS84_earth_major_equatorial_radius_of_ellipsoid * WGS84_earth_major_equatorial_radius_of_ellipsoid )
	one_f = ( 1 - WGS84_earth_ellipsoid_flattening )

	result = 0.0
	for i in range( 0 , integration_steps - 1 ):
		x = x1 + dx * i
		dydx_part1 = ( WGS84_earth_major_equatorial_radius_of_ellipsoid * one_f * math.sqrt( ( 1.0 - ( ( x + dx ) * ( x + dx ) ) / a_squared ) ) )
		dydx_part2 = ( WGS84_earth_major_equatorial_radius_of_ellipsoid * one_f * math.sqrt( ( 1.0 - ( x * x ) / a_squared ) ) )
		dydx = ( dydx_part1 - dydx_part2 ) / dx
		result += adx * math.sqrt( 1.0 + dydx * dydx )
	result_in_kilometers = result * 0.001
	return result_in_kilometers

def compute_distance_with_flattening( lat1 , lon1 , elevation_1_in_km , lat2 , lon2 , elevation_2_in_km ):
	earth_radius_in_km = 6371
	WGS84_earth_major_equatorial_radius_of_ellipsoid = 6378137
	WGS84_earth_ellipsoid_flattening = float( 1.0 / 298.257223563 )

	lat1_radians = lat1 * ( math.pi/180.0 )
	lon1_radians = lon1 * ( math.pi/180.0 )
	lat2_radians = lat2 * ( math.pi/180.0 )
	lon2_radians = lon2 * ( math.pi/180.0 )

	sin_of_lat1 = math.sin( lat1_radians )
	sin_of_lat2 = math.sin( lat2_radians )

	F = ( lat1_radians + lat2_radians ) / 2.0
	G = ( lat1_radians - lat2_radians ) / 2.0
	L = ( lon1_radians - lon2_radians ) / 2.0

	sing = math.sin( G )
	cosl = math.cos( L )
	cosf = math.cos( F )
	sinl = math.sin( L )
	sinf = math.sin( F )
	cosg = math.cos( G )

	S = sing*sing*cosl*cosl + cosf*cosf*sinl*sinl
	C = cosg*cosg*cosl*cosl + sinf*sinf*sinl*sinl
	W = math.atan2( math.sqrt( S ) , math.sqrt( C ) )
	R = math.sqrt( ( S * C ) ) / W
	H1 = ( 3 * R - 1.0 ) / ( 2.0 * C )
	H2 = ( 3 * R + 1.0 ) / ( 2.0 * S )
	D = 2 * W * earth_radius_in_km;
	return ( D * ( 1 + WGS84_earth_ellipsoid_flattening * H1 * sinf*sinf*cosg*cosg - WGS84_earth_ellipsoid_flattening*H2*cosf*cosf*sing*sing ) )

# https://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/Geographic-Distance-and-Azimuth-Calculations.htm
# Radius of Earth is in Kilometers, so Altitude Values must Match Units
def compute_distance_with_elevation( lat1 , lon1 , elevation_1_in_km , lat2 , lon2 , elevation_2_in_km ):
	earth_radius_in_km = 6371
	x1 = ( earth_radius_in_km + elevation_1_in_km ) * ( math.cos( lat1 ) * math.cos( lon1 ) )
	y1 = ( earth_radius_in_km + elevation_1_in_km ) * ( math.sin( lat1 ) )
	z1 = ( earth_radius_in_km + elevation_1_in_km ) * ( math.cos( lat1 ) * math.sin( lon1 ) )

	x2 = ( earth_radius_in_km + elevation_2_in_km ) * ( math.cos( lat2 ) * math.cos( lon2 ) )
	y2 = ( earth_radius_in_km + elevation_2_in_km ) * ( math.sin( lat2 ) )
	z2 = ( earth_radius_in_km + elevation_2_in_km ) * ( math.cos( lat2 ) * math.sin( lon2 ) )

	dx = ( x2 - x1 )
	dy = ( y2 - y1 )
	dz = ( z2 - z1 )
	distance = math.sqrt( dx*dx + dy*dy + dz*dz )
	return distance

# https://www.geeksforgeeks.org/haversine-formula-to-find-distance-between-two-points-on-a-sphere/
def compute_haversine_km_distance( lat1 , lon1 , lat2 , lon2 ):
	dLat = ( lat2 - lat1 ) * ( math.pi/180.0 )
	dLon = ( lon2 - lon1 ) * ( math.pi/180.0 )
	lat1 = (lat1) * ( math.pi/180.0 )
	lat2 = ( lat2 ) * ( math.pi/180.0 )
	angle = ( pow( math.sin( dLat/2 ) , 2 ) +  pow( math.sin( dLon/2 ) , 2 ) * math.cos( lat1 ) * math.cos( lat2 ) );
	earth_radius = 6371
	intermediate_c = ( 2 * math.asin( math.sqrt( angle ) ) )
	answer = ( earth_radius * intermediate_c )
	return answer

def location_info( location_1_name , location_2_name ):
	location_1 = geolocator.geocode( location_1_name )
	location_1_DMS = [ decimal_degrees_to_dms( location_1.latitude , 'lat' ) , decimal_degrees_to_dms( location_1.longitude , 'lon' ) ]
	location_2 = geolocator.geocode( location_2_name )
	location_2_DMS = [ decimal_degrees_to_dms( location_2.latitude , 'lat' ) , decimal_degrees_to_dms( location_2.longitude , 'lon' ) ]
	haversine_km_distance = compute_haversine_km_distance( location_1.latitude , location_1.longitude , location_2.latitude , location_2.longitude )
	great_circle_meter_distance = great_circle( ( location_1.latitude , location_1.longitude ) , ( location_2.latitude , location_2.longitude ) ).meters
	geodesic_meter_distance = geodesic( ( location_1.latitude , location_1.longitude ) , ( location_2.latitude , location_2.longitude ) ).meters
	vincenty_meter_distance = vincenty( ( location_1.latitude , location_1.longitude ) , ( location_2.latitude , location_2.longitude ) ).meters
	delta_altitudes = abs( location_1.altitude - location_2.altitude )
	distance_at_30000_feet = compute_distance_with_elevation( location_1.latitude , location_1.longitude , 9.144 , location_2.latitude , location_2.longitude , 9.144 )
	distance_with_flattening_in_km = compute_distance_with_flattening( location_1.latitude , location_1.longitude , 9.144 , location_2.latitude , location_2.longitude , 9.144 )
	print( f"{location_1_name} Latitude = {location_1.latitude}" )
	print( f"{location_1_name} Longitude = {location_1.longitude}" )
	print( f"{location_1_name} Latitude DMS = {location_1_DMS[0]}" )
	print( f"{location_1_name} Longitude DMS = {location_1_DMS[1]}" )
	print( f"{location_1_name} Altitude = {location_1.altitude}" )
	print( "" )
	print( f"{location_2_name} Latitude = {location_2.latitude}" )
	print( f"{location_2_name} Longitude = {location_2.longitude}" )
	print( f"{location_2_name} Latitude DMS = {location_2_DMS[0]}" )
	print( f"{location_2_name} Longitude DMS = {location_2_DMS[1]}" )
	print( f"{location_2_name} Altitude = {location_2.altitude}" )
	print( "" )
	print( f"Haversine KM Distance = {haversine_km_distance}" )
	print( f"Great Circle Meter Distance = {great_circle_meter_distance}" )
	print( f"Geodesic Meter Distance = {geodesic_meter_distance}" )
	print( f"Vincenty Meter Distance = {vincenty_meter_distance}" )
	print( f"KM Distance At 30,000 Feet = {distance_at_30000_feet}" )
	print( f"KM Distance With Flattening = {distance_with_flattening_in_km}" )
	return haversine_km_distance

#location_info( "Dayton Ohio" , "NorthPole" )
location_info( "Dayton Ohio Airport" , "San Francisco International Airport" )