Feat/geo

pyproject.toml

@@ -39,6 +39,10 @@ warn_return_any = true
 warn_unused_configs = true
 disallow_untyped_defs = true
 
+[tool.mypy.overrides]      
+module = "tests.*"
+disallow_untyped_defs = false
+
 [tool.ruff]
 line-length = 88
 target-version = "py311"

src/floatmatcher/geo.py

@@ -0,0 +1,20 @@
+# geo.py : all geography transformations
+
+# imports
+import numpy as np 
+from numpy.typing import ArrayLike, NDArray
+from .constants import EARTH_RADIUS_KM
+
+# functions 
+#_______
+
+def lonlat_to_xyz(lon: ArrayLike, lat: ArrayLike) -> NDArray[np.float64]:
+    """Convert lon/lat to xyz on the shpere"""
+    lon_rad = np.radians(lon)
+    lat_rad = np.radians(lat)
+    x = EARTH_RADIUS_KM * np.cos(lat_rad) * np.cos(lon_rad)
+    y = EARTH_RADIUS_KM * np.cos(lat_rad) * np.sin(lon_rad)
+    z = EARTH_RADIUS_KM * np.sin(lat_rad)
+    return np.column_stack([x, y, z])
+
+

tests/test_geo.py

@@ -0,0 +1,58 @@
+# tests/test_geo.py
+
+import numpy as np
+import numpy.testing as npt
+from floatmatcher.constants import EARTH_RADIUS_KM
+from floatmatcher.geo import (
+    lonlat_to_xyz
+)
+
+
+# ───────────── lonlat_to_xyz ─────────────
+
+def test_reference_points():
+    """ set special points to verify the computations """
+    lon = [0, 90, 0, 0]
+    lat = [0, 0, 90, -90]
+    xyz = lonlat_to_xyz(lon, lat)
+
+    expected = np.array([
+        [EARTH_RADIUS_KM, 0, 0],     # equator / Greenwich
+        [0, EARTH_RADIUS_KM, 0],     # equator / 90°E
+        [0, 0, EARTH_RADIUS_KM],     # north pole
+        [0, 0, -EARTH_RADIUS_KM],    # south pole
+    ])
+    npt.assert_allclose(xyz, expected, atol=1e-9)
+
+
+def test_output_shape():
+    """N points entry -> output (N, 3)"""
+    xyz = lonlat_to_xyz([0, 10, 20, 30], [0, 10, 20, 30])
+    assert xyz.shape == (4, 3)
+
+
+def test_points_lie_on_sphere():
+    """every point is at the same distance from the center"""
+    rng = np.random.default_rng(0)
+    lon = rng.uniform(-180, 180, size=100)
+    lat = rng.uniform(-90, 90, size=100)
+    xyz = lonlat_to_xyz(lon, lat)
+
+    norms = np.linalg.norm(xyz, axis=1) # compute norm of each xyz lon/lat vector 
+    npt.assert_allclose(norms, EARTH_RADIUS_KM, atol=1e-9)
+
+
+def test_antimeridian_points_are_close():
+    """2 points near -180 and 180 are close in 3D"""
+    xyz = lonlat_to_xyz([179, -179], [0, 0])
+    dist = np.linalg.norm(xyz[0] - xyz[1])
+    # equator: 2° ~222 km 
+    assert dist < 250
+
+
+def test_xyz_convention_invariant():
+    """each convention 0-360 and -180-180 must give the same xyz"""
+    a = lonlat_to_xyz([350.0], [30.0])   # 350 == -10 : same meridian
+    b = lonlat_to_xyz([-10.0], [30.0])
+    npt.assert_allclose(a, b, atol=1e-9)
+