Merge pull request #28 from joostlek/radius

Add methods to calculate bounding box from radius

Author: joostlek

src/python_opensky/opensky.py

@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import asyncio
+import math
 import socket
 from dataclasses import dataclass, field
 from datetime import datetime, timedelta, timezone
@@ -220,6 +221,130 @@ def _convert_state(state: list[Any]) -> dict[str, Any]:
 
         return dict(zip(keys, state, strict=True))
 
+    @staticmethod
+    # pylint: disable=too-many-locals
+    def calculate_point(
+        latitude: float,
+        longitude: float,
+        distance: float,
+        degrees: float,
+    ) -> tuple[float, float]:
+        """Calculate a point from an origin point, direction in degrees and distance."""
+        # ruff: noqa: N806
+        # pylint: disable=invalid-name
+        pi_d4 = math.atan(1.0)
+        two_pi = pi_d4 * 8.0
+        latitude = latitude * pi_d4 / 45.0
+        longitude = longitude * pi_d4 / 45.0
+        degrees = degrees * pi_d4 / 45.0
+        if degrees < 0.0:
+            degrees = degrees + two_pi
+        if degrees > two_pi:
+            degrees = degrees - two_pi
+        axis_a = 6378137
+        flattening = 1 / 298.257223563
+        axis_b = axis_a * (1.0 - flattening)
+        tan_u1 = (1 - flattening) * math.tan(latitude)
+        u1 = math.atan(tan_u1)
+        sigma1 = math.atan2(tan_u1, math.cos(degrees))
+        sinalpha = math.cos(u1) * math.sin(degrees)
+        cosalpha_sq = 1.0 - sinalpha * sinalpha
+        u2 = cosalpha_sq * (axis_a * axis_a - axis_b * axis_b) / (axis_b * axis_b)
+        A = 1.0 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
+        B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
+        # Starting with the approx
+        sigma = distance / (axis_b * A)
+        last_sigma = 2.0 * sigma + 2.0  # something impossible
+
+        # Iterate the following 3 eqs until no sig change in sigma
+        # two_sigma_m , delta_sigma
+        while abs((last_sigma - sigma) / sigma) > 1.0e-9:
+            two_sigma_m = 2 * sigma1 + sigma
+            delta_sigma = (
+                B
+                * math.sin(sigma)
+                * (
+                    math.cos(two_sigma_m)
+                    + (B / 4)
+                    * (
+                        math.cos(sigma)
+                        * (
+                            -1
+                            + 2 * math.pow(math.cos(two_sigma_m), 2)
+                            - (B / 6)
+                            * math.cos(two_sigma_m)
+                            * (-3 + 4 * math.pow(math.sin(sigma), 2))
+                            * (-3 + 4 * math.pow(math.cos(two_sigma_m), 2))
+                        )
+                    )
+                )
+            )
+            last_sigma = sigma
+            sigma = (distance / (axis_b * A)) + delta_sigma
+        phi2 = math.atan2(
+            (
+                math.sin(u1) * math.cos(sigma)
+                + math.cos(u1) * math.sin(sigma) * math.cos(degrees)
+            ),
+            (
+                (1 - flattening)
+                * math.sqrt(
+                    math.pow(sinalpha, 2)
+                    + pow(
+                        math.sin(u1) * math.sin(sigma)
+                        - math.cos(u1) * math.cos(sigma) * math.cos(degrees),
+                        2,
+                    ),
+                )
+            ),
+        )
+        lembda = math.atan2(
+            (math.sin(sigma) * math.sin(degrees)),
+            (
+                math.cos(u1) * math.cos(sigma)
+                - math.sin(u1) * math.sin(sigma) * math.cos(degrees)
+            ),
+        )
+        C = (flattening / 16) * cosalpha_sq * (4 + flattening * (4 - 3 * cosalpha_sq))
+        omega = lembda - (1 - C) * flattening * sinalpha * (
+            sigma
+            + C
+            * math.sin(sigma)
+            * (
+                math.cos(two_sigma_m)
+                + C * math.cos(sigma) * (-1 + 2 * math.pow(math.cos(two_sigma_m), 2))
+            )
+        )
+        lembda2 = longitude + omega
+        math.atan2(
+            sinalpha,
+            (
+                -math.sin(u1) * math.sin(sigma)
+                + math.cos(u1) * math.cos(sigma) * math.cos(degrees)
+            ),
+        )
+        phi2 = phi2 * 45.0 / pi_d4
+        lembda2 = lembda2 * 45.0 / pi_d4
+        return phi2, lembda2
+
+    @staticmethod
+    def get_bounding_box(
+        latitude: float,
+        longitude: float,
+        radius: float,
+    ) -> BoundingBox:
+        """Get bounding box from radius and a point."""
+        north = OpenSky.calculate_point(latitude, longitude, radius, 0)
+        east = OpenSky.calculate_point(latitude, longitude, radius, 90)
+        south = OpenSky.calculate_point(latitude, longitude, radius, 180)
+        west = OpenSky.calculate_point(latitude, longitude, radius, 270)
+        return BoundingBox(
+            min_latitude=min(north[0], south[0]) + latitude,
+            max_latitude=max(north[0], south[0]) + latitude,
+            min_longitude=min(east[1], west[1]) + longitude,
+            max_longitude=max(east[1], west[1]) + longitude,
+        )
+
     async def close(self) -> None:
         """Close open client session."""
         if self.session and self._close_session:

tests/test_radius.py

@@ -0,0 +1,23 @@
+"""Tests for the OpenSky Library."""
+
+
+from python_opensky import (
+    OpenSky,
+)
+
+
+async def test_calculating_bounding_box() -> None:
+    """Test calculating bounding box."""
+    bounding_box = OpenSky.get_bounding_box(0.0, 0.0, 25000)
+    assert bounding_box.min_latitude == -0.22609235747829648
+    assert bounding_box.max_latitude == 0.22609235747829648
+    assert bounding_box.min_longitude == -0.22457882102988042
+    assert bounding_box.max_longitude == 0.22457882102988042
+
+
+async def test_calculating_direction() -> None:
+    """Test calculating direction."""
+    second_point = OpenSky.calculate_point(0.0, 0.0, 25000.0, -180)
+    assert second_point == (-0.22609235747829648, 2.7503115231199028e-17)
+    second_point = OpenSky.calculate_point(0.0, 0.0, 25000.0, 361)
+    assert second_point == (0.22605792234324162, 0.003919461063277522)