Repository: ecmwf/earthkit-data
Pull request #229 (2023-10-13 17:25:39)

WIP: Implement haversine nearest point method
First implementation.

See the documentation in the `geo/distance.py` code.

This is a working example to extract time series for a given point:

```python
from earthkit.data import from_source
from earthkit.data.geo import nearest_point_haversine

ds = from_source("file", "tests/data/t_time_series.grib").sel(param="t")

latlon = ds.to_latlon()
lat = latlon["lat"]
lon = latlon["lon"]

p_ref = (51.45, -0.97)
idx = nearest_point_haversine(p_ref, (lat, lon))
```
12

```python
ds.values[:,idx]
```
array([280.44058228, 280.31297302, 280.2789917 , 280.08499146,
       276.38993835])



Patch:
diff --git a/earthkit/data/core/constants.py b/earthkit/data/core/constants.py
index ff77166..d62f961 100644
--- a/earthkit/data/core/constants.py
+++ b/earthkit/data/core/constants.py
@@ -8,3 +8,21 @@
 #
 
 DATETIME = "valid_datetime"
+
+
+"""
+Collection of constants in SI units.
+"""
+
+R_earth = 6371229
+r"""Average radius of the Earth [:math:`m`]. See [IFS-CY47R3-PhysicalProcesses]_
+ (Chapter 12)."""
+
+full_circle = 360
+r"""Full circle in degrees"""
+
+north = 90
+r"""Latitude of the north pole in degrees"""
+
+south = -90
+r"""Latitude of the south pole in degrees"""
diff --git a/earthkit/data/geo/__init__.py b/earthkit/data/geo/__init__.py
new file mode 100644
index 0000000..9e5c263
--- /dev/null
+++ b/earthkit/data/geo/__init__.py
@@ -0,0 +1,10 @@
+# (C) Copyright 2020 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+
+from .distance import haversine_distance, nearest_point_haversine  # noqa
diff --git a/earthkit/data/geo/distance.py b/earthkit/data/geo/distance.py
new file mode 100644
index 0000000..6f4564c
--- /dev/null
+++ b/earthkit/data/geo/distance.py
@@ -0,0 +1,125 @@
+# (C) Copyright 2020 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+
+import numpy as np
+
+from earthkit.data.core import constants
+
+
+def regulate_lat(lat):
+    return np.where(np.abs(lat) > constants.north, np.nan, lat)
+
+
+def haversine_distance(p1, p2):
+    """Compute haversine distance between two (sets of) points on Earth.
+
+    Parameters
+    ----------
+    p1: pair of array-like
+        Locations of the first points. The first item specifies the latitudes,
+        the second the longitudes (degrees)
+    p2: pair of array-like
+        Locations of the second points. The first item specifies the latitudes,
+        the second the longitudes (degrees)
+
+    Returns
+    -------
+    number or ndarray
+        Spherical distance on the surface in Earth (m)
+
+    Either ``p1`` or ``p2`` must be a single point.
+
+    Examples
+    --------
+    Compute the distance between Reading and Bologna.
+    >>> from earthkit.data.geo import haversine_distance
+    >>> p1 = (51.45, -0.97)
+    >>> p2 = (44.49, 11.34)
+    >>> haversine_distance(p1, p2)
+    1196782.5785709629
+
+    Compute the distance from Reading to Bologna and Bonn.
+    >>> from earthkit.data.geo import haversine_distance
+    >>> p1 = (51.45, -0.97)
+    >>> p_lat = [44.49, 50.73]
+    >>> p_lon = [11.34, 7.90]
+    >>> haversine_distance(p1, (p_lat, p_lon))
+    array([1196782.57857096,  624273.19519049])
+
+    """
+    lat1 = np.asarray(p1[0])
+    lon1 = np.asarray(p1[1])
+    lat2 = np.asarray(p2[0])
+    lon2 = np.asarray(p2[1])
+
+    if lat1.shape != lon1.shape:
+        raise ValueError(
+            f"haversine_distance: lat and lon in p1 must have the same shape! {lat1.shape} != {lon1.shape}"
+        )
+
+    if lat2.shape != lon2.shape:
+        raise ValueError(
+            f"haversine_distance: lat and lon in p2 must have the same shape! {lat2.shape} != {lon2.shape}"
+        )
+
+    if lat1.size != 1 and lat2.size != 1:
+        raise ValueError("haversine_distance: either p1 or p2 must be a single point")
+
+    lat1 = regulate_lat(lat1)
+    lat2 = regulate_lat(lat2)
+
+    lat1, lon1, lat2, lon2 = map(np.deg2rad, [lat1, lon1, lat2, lon2])
+    d_lon = lon2 - lon1
+    d_lat = lat2 - lat1
+
+    a = np.sqrt(
+        np.sin(d_lat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(d_lon / 2) ** 2
+    )
+    c = 2 * np.arcsin(a)
+    distance = constants.R_earth * c
+
+    return distance
+
+
+def nearest_point_haversine(ref_point, points):
+    """Find the index of the nearest point to ``ref_point`` in a set of ``points`` using the
+       haversine distance formula.
+
+    Parameters
+    ----------
+    point_ref: pair of numbers
+        Latitude and longitude coordinate of the reference point (degrees)
+    points: pair of array-like
+        Locations of the set of points from which the nearest to
+        ``point_ref`` is to be found. The first item specifies the latitudes,
+        the second the longitudes (degrees)
+
+    Returns
+    -------
+    number
+        Index of the nearest point from ``points`` to ``ref_point``
+
+    Examples
+    --------
+    >>> from earthkit.data.geo import nearest_point_haversine
+    >>> p_ref = (51.45, -0.97)
+    >>> p_lat = [44.49, 50.73, 50.1]
+    >>> p_lon = [11.34, 7.90, -8.1]
+    >>> nearest_point_haversine(p_ref, (p_lat, p_lon))
+    2
+
+    """
+    if np.asarray(ref_point[0]).size != 1 or np.asarray(ref_point[1]).size != 1:
+        raise ValueError("nearest_point_haversine: ref_point must be a single point")
+
+    distance = haversine_distance(ref_point, points)
+    index = np.nanargmin(distance)
+    if isinstance(index, np.ndarray):
+        return index[0]
+    return index


Test patch:
diff --git a/tests/geo/test_geo.py b/tests/geo/test_geo.py
new file mode 100644
index 0000000..99a42e5
--- /dev/null
+++ b/tests/geo/test_geo.py
@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+
+# (C) Copyright 2020 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+
+
+import numpy as np
+import pytest
+
+from earthkit.data.geo import haversine_distance, nearest_point_haversine
+
+
+def test_haversine_distance_single_point():
+    p_ref = (48.0, 20)
+    d_ref = 5675597.92279149
+    p = (0, 0)
+    d = haversine_distance(p_ref, p)
+    assert np.isclose(d, d_ref)
+
+
+@pytest.mark.parametrize("dlon_ref", [-720, -360, 0, 360, 720])
+@pytest.mark.parametrize("dlon_p", [-720, -360, 0, 360, 720])
+@pytest.mark.parametrize(
+    "p_ref,d_ref",
+    [
+        (
+            (0, 0),
+            [
+                0.0,
+                10007903.11036912,
+                10007903.11036912,
+                20015806.22073824,
+                10007903.11036912,
+                10007903.11036912,
+                5675597.92279149,
+                5675597.92279149,
+                5675597.92279149,
+                5675597.92279149,
+                np.nan,
+            ],
+        ),
+        (
+            (-90, 0),
+            [
+                10007903.11036912,
+                10007903.11036912,
+                10007903.11036912,
+                10007903.11036912,
+                20015806.22073824,
+                0.0,
+                15345451.43589932,
+                15345451.43589932,
+                4670354.78483892,
+                4670354.78483892,
+                np.nan,
+            ],
+        ),
+        (
+            (48.0, 20),
+            [
+                5675597.92279149,
+                8536770.52794796,
+                11479035.69279028,
+                14340208.29794676,
+                4670354.78483892,
+                15345451.43589932,
+                0.0,
+                2942265.16484232,
+                11351195.84558298,
+                10675096.6510604,
+                np.nan,
+            ],
+        ),
+    ],
+)
+def test_haversine_distance(dlon_ref, dlon_p, p_ref, d_ref):
+    lats = np.array([0.0, 0, 0, 0, 90, -90, 48, 48, -48, -48, np.nan])
+    lons = np.array([0, 90, -90, 180, 0, 0, 20, -20, -20, 20, 1.0]) + dlon_p
+    p_ref = (p_ref[0], p_ref[1] + dlon_ref)
+    d_ref = np.array(d_ref)
+
+    d = haversine_distance(p_ref, (lats, lons))
+    assert np.allclose(
+        d[:-1], d_ref[:-1]
+    ), f"p_ref={p_ref},dlon_ref={dlon_ref},dlon_p={dlon_p}"
+
+
+# lat is out of range
+@pytest.mark.parametrize("p_ref", [((90.00001, 0)), (((-90.00001, 0)))])
+def test_haversine_distance_invalid(p_ref):
+    p = (-48, 12)
+    d = haversine_distance(p_ref, p)
+    assert np.isnan(d)
+
+
+@pytest.mark.parametrize(
+    "p_ref,index_ref",
+    [
+        ((0, 0), 0),
+        ((15, 22), 0),
+        ((44, 10), 6),
+        ((44, -10), 7),
+        ((89, -120), 4),
+        ((-50, -18), 8),
+        ((-50, 18), 9),
+    ],
+)
+def test_haversine_nearest(p_ref, index_ref):
+    lats = np.array([0.0, 0, 0, 0, 90, -90, 48, 48, -48, -48, np.nan])
+    lons = np.array([0, 90, -90, 180, 0, 0, 20, -20, -20, 20, 1.0])
+
+    index = nearest_point_haversine(p_ref, (lats, lons))
+    assert index == index_ref, p_ref
+
+
+def test_haversine_nearest_invalid():
+    # p_ref must be a single point
+    p_ref = ([15, 21], [22, 7])
+    p = ([0.0, 0], [0, 90])
+
+    with pytest.raises(ValueError):
+        nearest_point_haversine(p_ref, p)
+
+
+if __name__ == "__main__":
+    from earthkit.data.testing import main
+
+    main(__file__)