diff --git a/ogn/commands/dbutils.py b/ogn/commands/dbutils.py index 070d7fc..99bfaf3 100644 --- a/ogn/commands/dbutils.py +++ b/ogn/commands/dbutils.py @@ -6,7 +6,7 @@ from sqlalchemy.sql import null from sqlalchemy.orm import sessionmaker from ogn.model import AircraftBeacon, ReceiverBeacon -from ogn.utils import wgs84_to_sphere +from ogn.utils import haversine_distance os.environ.setdefault('OGN_CONFIG_MODULE', 'config.default') @@ -33,9 +33,16 @@ def update_receiver_childs(name): AircraftBeacon.radius == null())) for aircraft_beacon in aircraft_beacons_query.all(): - [radius, theta, phi] = wgs84_to_sphere(last_receiver_beacon, - aircraft_beacon) - aircraft_beacon.radius = radius + location0 = (last_receiver_beacon.latitude, last_receiver_beacon.longitude) + location1 = (aircraft_beacon.latitude, aircraft_beacon.longitude) + alt0 = last_receiver_beacon.altitude + alt1 = aircraft_beacon.altitude + + (flat_distance, phi) = haversine_distance(location0, location1) + theta = atan2(alt1 - alt0, flat_distance) * 180 / pi + distance = sqrt(flat_distance**2 + (alt1 - alt0)**2) + + aircraft_beacon.radius = distance aircraft_beacon.theta = theta aircraft_beacon.phi = phi session.commit() diff --git a/ogn/utils.py b/ogn/utils.py index 698ef68..a8f546d 100644 --- a/ogn/utils.py +++ b/ogn/utils.py @@ -62,22 +62,16 @@ def get_country_code(latitude, longitude): return country_code -def wgs84_to_sphere(receiver_beacon, aircraft_beacon): - from math import pi, asin, sqrt, sin, cos, atan2 - deg2rad = pi / 180 - rad2deg = 180 / pi +def haversine_distance(location0, location1): + from math import asin, sqrt, sin, cos, atan2, radians, degrees - lat1 = receiver_beacon.latitude * deg2rad - lon1 = receiver_beacon.longitude * deg2rad - alt1 = receiver_beacon.altitude + lat0 = radians(location0[0]) + lon0 = radians(location0[1]) + lat1 = radians(location1[0]) + lon1 = radians(location1[1]) - lat2 = aircraft_beacon.latitude * deg2rad - lon2 = aircraft_beacon.longitude * deg2rad - alt2 = aircraft_beacon.altitude + distance = 6366000 * 2 * asin(sqrt((sin((lat0 - lat1) / 2))**2 + cos(lat0) * cos(lat1) * (sin((lon0 - lon1) / 2))**2)) + print(distance) + phi = degrees(atan2(sin(lon0 - lon1) * cos(lat1), cos(lat0) * sin(lat1) - sin(lat0) * cos(lat1) * cos(lon0 - lon1))) - distance = 6366000 * 2 * asin(sqrt((sin((lat1 - lat2) / 2))**2 + cos(lat1) * cos(lat2) * (sin((lon1 - lon2) / 2))**2)) - theta = atan2(alt2 - alt1, distance) * rad2deg - phi = atan2(sin(lon1 - lon2) * cos(lat2), cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon1 - lon2)) * rad2deg - - radius = sqrt(distance**2 + (alt2 - alt1)**2) - return radius, theta, phi + return distance, phi diff --git a/tests/test_utils.py b/tests/test_utils.py index 8a18173..c722b67 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,7 +1,7 @@ import unittest -from ogn.utils import get_ddb, get_trackable, get_country_code, wgs84_to_sphere -from ogn.model import AddressOrigin, Beacon +from ogn.utils import get_ddb, get_trackable, get_country_code, haversine_distance +from ogn.model import AddressOrigin class TestStringMethods(unittest.TestCase): @@ -44,51 +44,27 @@ class TestStringMethods(unittest.TestCase): country_code = get_country_code(latitude, longitude) self.assertEqual(country_code, None) - def test_wgs84_to_sphere(self): - receiver_beacon = Beacon() - receiver_beacon.latitude = 0 - receiver_beacon.longitude = 0 - receiver_beacon.altitude = 0 - + def test_haversine_distance(self): # delta: one latitude degree - aircraft_beacon = Beacon() - aircraft_beacon.latitude = -1 - aircraft_beacon.longitude = 0 - aircraft_beacon.altitude = 0 - [radius, theta, phi] = wgs84_to_sphere(receiver_beacon, aircraft_beacon) - self.assertAlmostEqual(radius, 60 * 1852, -2) - self.assertEqual(theta, 0) + location0 = (0, 0) + location1 = (-1, 0) + + (distance, phi) = haversine_distance(location0, location1) + self.assertAlmostEqual(distance, 60 * 1852, -2) self.assertEqual(phi, 180) # delta: one longitude degree at the equator - aircraft_beacon.latitude = 0 - aircraft_beacon.longitude = -1 - aircraft_beacon.altitude = 0 - [radius, theta, phi] = wgs84_to_sphere(receiver_beacon, aircraft_beacon) - self.assertAlmostEqual(radius, 60 * 1852, -2) - self.assertEqual(theta, 0) + location0 = (0, 0) + location1 = (0, -1) + + (distance, phi) = haversine_distance(location0, location1) + self.assertAlmostEqual(distance, 60 * 1852, -2) self.assertEqual(phi, 90) - # delta: 1000m altitude - aircraft_beacon.latitude = 0 - aircraft_beacon.longitude = 0 - aircraft_beacon.altitude = 1000 - [radius, theta, phi] = wgs84_to_sphere(receiver_beacon, aircraft_beacon) - self.assertAlmostEqual(radius, 1000, 3) - self.assertEqual(theta, 90) - self.assertEqual(phi, 0) + # delta: 29000m + location0 = (48.865, 9.2225) + location1 = (48.74435, 9.578) - # receiver - receiver_beacon.latitude = 48.865 - receiver_beacon.longitude = 9.2225 - receiver_beacon.altitude = 574 - - # aircraft beacon - aircraft_beacon.latitude = 48.74435 - aircraft_beacon.longitude = 9.578 - aircraft_beacon.altitude = 929 - - [radius, theta, phi] = wgs84_to_sphere(receiver_beacon, aircraft_beacon) - self.assertAlmostEqual(radius, 29265.6035812215, -1) - self.assertAlmostEqual(theta, 0.694979846308314, 5) + (distance, phi) = haversine_distance(location0, location1) + self.assertAlmostEqual(distance, 29265.6035812215, -1) self.assertAlmostEqual(phi, -117.1275408121, 5)