Add coordinate and wind rotation (#10)

* Add point and wind rotation

docs/conf.py

@@ -29,6 +29,7 @@
 # ones.
 extensions = [
     "sphinx_rtd_theme",
+    "nbsphinx",
     "sphinx.ext.autodoc",
     "sphinx.ext.napoleon",
     "autoapi.extension",
@@ -45,7 +46,6 @@
     "undoc-members",
     "show-inheritance",
     "show-module-summary",
-    "imported-members",
     "inherited-members",
 ]
 autoapi_root = "_api"

docs/examples/index.rst

@@ -0,0 +1,11 @@
+.. _examples:
+
+Examples
+============
+
+Here is a list of example notebooks to illustrate how to use earthkit-geo.
+
+.. toctree::
+    :maxdepth: 1
+
+    rotate.ipynb

docs/index.rst

@@ -13,9 +13,9 @@ Welcome to earthkit-geo's documentation
 .. toctree::
    :maxdepth: 1
    :caption: Examples
-   :titlesonly:
 
-   examples
+   examples/index
+
 
 .. toctree::
    :maxdepth: 1

docs/references.rst

@@ -4,3 +4,19 @@ References
 .. [IFS-CY47R3-PhysicalProcesses]
 
 IFS Documentation CY47R3 - Part IV Physical processes, (2021). URL: https://www.ecmwf.int/en/elibrary/20198-ifs-documentation-cy47r3-part-iv-physical-processes
+
+
+
+.. [ECEF]
+
+Earth-centred, Earth-fixed coordinate system. URL: https://en.wikipedia.org/wiki/Earth-centered,_Earth-fixed_coordinate_system
+
+
+.. [From_ECEF_to_geodetic_coordinates]
+
+URL: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_geodetic_coordinates
+
+
+.. [From_geodetic_to_ECEF_coordinates]
+
+URL: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

docs/requirements.txt

@@ -2,7 +2,8 @@
 ipykernel
 docutils
 Pygments>=2.6.1
-Sphinx
+Sphinx>=7.3.7
 sphinx-rtd-theme
+nbsphinx
 setuptools
 sphinx-autoapi

environment.yml

@@ -3,6 +3,7 @@ channels:
 - conda-forge
 - nodefaults
 dependencies:
+- pyproj
 - scipy
 - make
 - mypy
@@ -15,3 +16,4 @@ dependencies:
   - sphinx-autoapi>=3.0.0
 - sphinx_rtd_theme
 - sphinxcontrib-apidoc
+- nbsphinx

pyproject.toml

@@ -19,6 +19,7 @@ classifiers = [
   "Operating System :: OS Independent"
 ]
 dependencies = [
+  "pyproj",
   "scipy"
 ]
 description = "Geospatial computations"
@@ -37,7 +38,7 @@ test = [
 [project.urls]
 Documentation = "https://earthkit-geo.readthedocs.io/"
 Homepage = "https://github.com/ecmwf/earthkit-geo/"
-Issues = "https://github.com/ecmwf/earthkit-geo.issues"
+Issues = "https://github.com/ecmwf/earthkit-geo/issues"
 Repository = "https://github.com/ecmwf/earthkit-geo/"
 
 [tool.coverage.run]

src/earthkit/geo/constants.py

@@ -11,5 +11,14 @@
 Collection of constants in SI units.
 """
 
-north = 90
+NORTH_POLE_LAT = 90
 r"""Latitude of the north pole in degrees"""
+
+SOUTH_POLE_LAT = -90
+r"""Latitude of the south pole in degrees"""
+
+FULL_ANGLE = 360.0
+r"""Full angle in degrees"""
+
+STRAIGHT_ANGLE = 180.0
+r"""Half angle in degrees"""

src/earthkit/geo/coord.py

@@ -0,0 +1,92 @@
+# (C) Copyright 2024 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 . import constants
+
+
+def _normalise_lon(lon, minimum):
+    while lon < minimum:
+        lon += constants.FULL_ANGLE
+
+    while lon >= minimum + constants.FULL_ANGLE:
+        lon -= constants.FULL_ANGLE
+
+    return lon
+
+
+def xyz_to_latlon(x, y, z):
+    """Convert from earth-centred, earth-fixed ([ECEF]_) to geodetic coordinates.
+    See [From_ECEF_to_geodetic_coordinates]_.
+
+    Parameters
+    ----------
+    x: float or ndarray
+        x-component of the ECEF coordinates.
+    y: float or ndarray
+        y-component of the ECEF coordinates.
+    z: float or ndarray
+        z-component of the ECEF coordinates.
+
+    Returns
+    -------
+    ndarray
+        Latitude (degrees).
+    ndarray
+        Longitude (degrees).
+
+
+    It is assumed that the Earth is a sphere of radius 1.
+    """
+    return (
+        np.rad2deg(np.arcsin(np.minimum(1.0, np.maximum(-1.0, z)))),
+        np.rad2deg(np.arctan2(y, x)),
+    )
+
+
+def latlon_to_xyz(lat, lon):
+    """Convert from geodetic to earth-centred, earth-fixed ([ECEF]_) coordinates.
+    See [From_geodetic_to_ECEF_coordinates]_.
+
+    Parameters
+    ----------
+    lat: float or ndarray
+        Latitude (degrees).
+    lon: float or ndarray
+        Longitude (degrees).
+
+    Returns
+    -------
+    ndarray
+        x-component of the ECEF coordinates.
+    ndarray
+        y-component of the ECEF coordinates.
+    ndarray
+        z-component of the ECEF coordinates.
+
+
+    It is assumed that the Earth is a sphere of radius 1. It is also assumed the
+    geodetic coordinate h = 0.
+    """
+    lat = np.asarray(lat)
+    lon = np.asarray(lon)
+    phi = np.deg2rad(lat)
+    lda = np.deg2rad(lon)
+
+    cos_phi = np.cos(phi)
+    cos_lda = np.cos(lda)
+    sin_phi = np.sin(phi)
+    sin_lda = np.sin(lda)
+
+    x = cos_phi * cos_lda
+    y = cos_phi * sin_lda
+    z = sin_phi
+
+    return x, y, z

src/earthkit/geo/distance.py

@@ -10,11 +10,12 @@
 import numpy as np
 
 from . import constants
+from .coord import latlon_to_xyz
 from .figure import IFS_SPHERE, UNIT_SPHERE
 
 
 def _regulate_lat(lat):
-    return np.where(np.abs(lat) > constants.north, np.nan, lat)
+    return np.where(np.abs(lat) > constants.NORTH_POLE_LAT, np.nan, lat)
 
 
 def haversine_distance(p1, p2, figure=IFS_SPHERE):
@@ -152,19 +153,6 @@ def nearest_point_haversine(ref_points, points, figure=IFS_SPHERE):
     return (np.array(res_index), figure.scale(np.array(res_distance)))
 
 
-def _latlon_to_xyz(lat, lon):
-    """Works on the unit sphere."""
-    lat = np.asarray(lat)
-    lon = np.asarray(lon)
-    lat = np.radians(lat)
-    lon = np.radians(lon)
-    x = np.cos(lat) * np.cos(lon)
-    y = np.cos(lat) * np.sin(lon)
-    z = np.sin(lat)
-
-    return x, y, z
-
-
 def _chordlength_to_arclength(chord_length):
     """
     Convert 3D (Euclidean) distance to great circle arc length
@@ -199,7 +187,7 @@ def __init__(self, lats, lons):
             lats = lats[mask]
             lons = lons[mask]
 
-        x, y, z = _latlon_to_xyz(lats, lons)
+        x, y, z = latlon_to_xyz(lats, lons)
         v = np.column_stack((x, y, z))
         self.tree = KDTree(v)
 
@@ -232,7 +220,7 @@ def nearest_point(self, ref_points, figure=IFS_SPHERE):
 
         """
         lat, lon = ref_points
-        x, y, z = _latlon_to_xyz(lat, lon)
+        x, y, z = latlon_to_xyz(lat, lon)
         points = np.column_stack((x, y, z))
 
         # find the nearest point