Add timezone and observing night calculations

pyproject.toml

@@ -45,6 +45,7 @@ dependencies = [
   "naif-eop-predict",
   "naif-eop-historical",
   "naif-earth-itrf93",
+  "timezonefinder",
 ]
 
 [build-system]

src/adam_core/constants.py

@@ -14,11 +14,19 @@
 
 
 class _Constants:
-    def __init__(self, C=None, MU=None, R_Earth=None, Obliquity=None):
+    def __init__(
+        self,
+        C=None,
+        MU=None,
+        R_EARTH_EQUATORIAL=None,
+        R_EARTH_POLAR=None,
+        OBLIQUITY=None,
+    ):
         self.C = C
         self.MU = MU
-        self.R_EARTH = R_Earth
-        self.OBLIQUITY = Obliquity
+        self.R_EARTH_EQUATORIAL = R_EARTH_EQUATORIAL
+        self.R_EARTH_POLAR = R_EARTH_POLAR
+        self.OBLIQUITY = OBLIQUITY
 
         # Transformation matrix from Equatorial J2000 to Ecliptic J2000
         self.TRANSFORM_EQ2EC = np.array(
@@ -40,8 +48,10 @@ def __init__(self, C=None, MU=None, R_Earth=None, Obliquity=None):
     # Standard Gravitational Parameter -- Sun :  au**3 / d**2 (0.29591220828411956E-03 -- DE441/DE440)
     "MU": 0.29591220828411956e-03,
     # Earth Equatorial Radius: au (6378.1363 km -- DE431/DE430)
-    "R_Earth": 6378.1363 / KM_P_AU,
+    "R_EARTH_EQUATORIAL": 6378.1363 / KM_P_AU,
+    # Earth Polar Radius: au (6356.7523 km)
+    "R_EARTH_POLAR": 6356.7523 / KM_P_AU,
     # Mean Obliquity at J2000: radians (84381.448 arcseconds -- DE431/DE430)
-    "Obliquity": 84381.448 * np.pi / (180.0 * 3600.0),
+    "OBLIQUITY": 84381.448 * np.pi / (180.0 * 3600.0),
 }
 DE44X = Constants = _Constants(**DE44X_CONSTANTS)

src/adam_core/observers/__init__.py

@@ -1,5 +1,6 @@
 # flake8: noqa: F401
 from .observers import Observers
 from .state import get_observer_state
+from .utils import calculate_observing_night
 
 __all__ = ["Observers", "get_observer_state"]

src/adam_core/observers/observers.py

@@ -8,20 +8,80 @@
 import pyarrow.compute as pc
 import quivr as qv
 from mpc_obscodes import mpc_obscodes
+from timezonefinder import TimezoneFinder
 from typing_extensions import Self
 
+from ..constants import Constants as c
 from ..coordinates.cartesian import CartesianCoordinates
 from ..coordinates.origin import OriginCodes
 from ..time import Timestamp
 
+R_EARTH_EQUATORIAL = c.R_EARTH_EQUATORIAL
+R_EARTH_POLAR = c.R_EARTH_POLAR
+E_EARTH = np.sqrt(1 - (R_EARTH_POLAR / R_EARTH_EQUATORIAL) ** 2)
 
-class ObservatoryGeodetics(qv.Table):
+
+class ObservatoryParallaxCoefficients(qv.Table):
     code = qv.LargeStringColumn()
     longitude = qv.Float64Column()
     cos_phi = qv.Float64Column()
     sin_phi = qv.Float64Column()
     name = qv.LargeStringColumn()
 
+    def lon_lat(self) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+        """
+        Return the longitude and latitude of the observatories in degrees.
+
+        This is only valid for Earth-based observatories.
+
+        Returns
+        -------
+        longitude : np.ndarray
+            The longitude of the observatories in degrees. In the range -180 to 180 degrees,
+            with positive values east of the prime meridian.
+        latitude : np.ndarray
+            The latitude of the observatories in degrees. In the range -90 to 90 degrees,
+            with positive values north of the equator.
+        """
+        # Filter out Space-based observatories
+        mask = pc.is_nan(self.longitude).to_numpy(zero_copy_only=False)
+
+        longitude = np.where(
+            mask, np.nan, self.longitude.to_numpy(zero_copy_only=False)
+        )
+        tan_phi_geo = np.where(
+            mask,
+            np.nan,
+            self.sin_phi.to_numpy(zero_copy_only=False)
+            / self.cos_phi.to_numpy(zero_copy_only=False),
+        )
+        latitude_geodetic = np.arctan(tan_phi_geo / (1 - E_EARTH**2))
+
+        # Scale longitude to -180 to 180
+        longitude = np.where(longitude > 180, longitude - 360, longitude)
+
+        return longitude, np.degrees(latitude_geodetic)
+
+    def timezone(self) -> npt.NDArray[np.str_]:
+        """
+        Return the timezone of the observatories in hours.
+
+        Returns
+        -------
+        timezone : np.ndarray
+            The timezone of the observatories in hours.
+        """
+        tf = TimezoneFinder()
+        lon, lat = self.lon_lat()
+        time_zones = np.array(
+            [
+                tz if not np.isnan(lon_i) else "None"
+                for lon_i, lat_i in zip(lon, lat)
+                for tz in [tf.timezone_at(lng=lon_i, lat=lat_i)]
+            ]
+        )
+        return time_zones
+
 
 # Read MPC extended observatory codes file
 # Ignore warning about pandas deprecation
@@ -39,7 +99,7 @@ class ObservatoryGeodetics(qv.Table):
     )
     OBSCODES.reset_index(inplace=True, names=["code"])
 
-OBSERVATORY_GEODETICS = ObservatoryGeodetics.from_kwargs(
+OBSERVATORY_PARALLAX_COEFFICIENTS = ObservatoryParallaxCoefficients.from_kwargs(
     code=OBSCODES["code"].values,
     longitude=OBSCODES["Longitude"].values,
     cos_phi=OBSCODES["cos"].values,
@@ -48,7 +108,7 @@ class ObservatoryGeodetics(qv.Table):
 )
 
 OBSERVATORY_CODES = {
-    x for x in OBSERVATORY_GEODETICS.code.to_numpy(zero_copy_only=False)
+    x for x in OBSERVATORY_PARALLAX_COEFFICIENTS.code.to_numpy(zero_copy_only=False)
 }
 
 

src/adam_core/observers/state.py

@@ -9,9 +9,9 @@
 from ..coordinates.origin import Origin, OriginCodes
 from ..time import Timestamp
 from ..utils.spice import get_perturber_state, setup_SPICE
-from .observers import OBSERVATORY_CODES, OBSERVATORY_GEODETICS
+from .observers import OBSERVATORY_CODES, OBSERVATORY_PARALLAX_COEFFICIENTS
 
-R_EARTH = c.R_EARTH
+R_EARTH_EQUATORIAL = c.R_EARTH_EQUATORIAL
 OMEGA_EARTH = 2 * np.pi / 0.997269675925926
 Z_AXIS = np.array([0, 0, 1])
 
@@ -77,17 +77,17 @@ def get_observer_state(
     # If it not an origin code then lets get the state vector
     else:
         # Get observatory geodetic information
-        geodetics = OBSERVATORY_GEODETICS.select("code", code)
-        if len(geodetics) == 0:
+        parallax_coeffs = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", code)
+        if len(parallax_coeffs) == 0:
             raise ValueError(
-                f"Observatory code '{code}' not found in the observatory geodetics table."
+                f"Observatory code '{code}' not found in the observatory parallax coefficients table."
             )
-        geodetics = geodetics.table.to_pylist()[0]
+        parallax_coeffs = parallax_coeffs.table.to_pylist()[0]
 
         # Unpack geodetic information
-        longitude = geodetics["longitude"]
-        cos_phi = geodetics["cos_phi"]
-        sin_phi = geodetics["sin_phi"]
+        longitude = parallax_coeffs["longitude"]
+        cos_phi = parallax_coeffs["cos_phi"]
+        sin_phi = parallax_coeffs["sin_phi"]
 
         if np.any(np.isnan([longitude, cos_phi, sin_phi])):
             err = (
@@ -118,7 +118,7 @@ def get_observer_state(
             )
 
             # Multiply pointing vector with Earth radius to get actual vector
-            o_vec_ITRF93 = np.dot(R_EARTH, o_hat_ITRF93)
+            o_vec_ITRF93 = np.dot(R_EARTH_EQUATORIAL, o_hat_ITRF93)
 
             # Warning! Converting times to ET will incur a loss of precision.
             epochs_et = times.rescale("tdb").et()
@@ -136,7 +136,6 @@ def get_observer_state(
                 # Nutation:  1980 IAU model, with IERS corrections due to Herring et al.
                 # True sidereal time using accurate values of TAI-UT1
                 # Polar motion
-
                 rotation_matrix = sp.pxform("ITRF93", frame_spice, epoch.as_py())
 
                 # Find indices of epochs that match the current unique epoch
@@ -148,7 +147,7 @@ def get_observer_state(
                 # Calculate the velocity (thank you numpy broadcasting)
                 rotation_direction = np.cross(o_hat_ITRF93, Z_AXIS)
                 v_obs[mask] = v_geo[mask] + rotation_matrix @ (
-                    -OMEGA_EARTH * R_EARTH * rotation_direction
+                    -OMEGA_EARTH * R_EARTH_EQUATORIAL * rotation_direction
                 )
 
         return CartesianCoordinates.from_kwargs(

src/adam_core/observers/tests/test_observers.py

@@ -1,9 +1,10 @@
+import numpy as np
 import pyarrow as pa
 import pyarrow.compute as pc
 import pytest
 
 from ...time import Timestamp
-from ..observers import Observers
+from ..observers import OBSERVATORY_PARALLAX_COEFFICIENTS, Observers
 
 
 @pytest.fixture
@@ -59,3 +60,53 @@ def test_Observers_from_codes_raises(codes_times) -> None:
         Observers.from_codes(codes[:3], times)
     with pytest.raises(ValueError, match="codes and times must have the same length."):
         Observers.from_codes(codes, times[:3])
+
+
+def test_ObservatoryParallaxCoefficients_lon_lat() -> None:
+    # Data taken from: https://en.wikipedia.org/wiki/List_of_observatory_codes
+    # From: https://geohack.toolforge.org/geohack.php?pagename=Zwicky_Transient_Facility&params=33.35731_N_116.85981_W_
+    I41 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "I41")
+    lon, lat = I41.lon_lat()
+    np.testing.assert_allclose(lon[0], -116.85981, atol=1e-4, rtol=0)
+    np.testing.assert_allclose(lat[0], 33.35731, atol=1e-4, rtol=0)
+
+    # From: https://geohack.toolforge.org/geohack.php?pagename=Vera_C._Rubin_Observatory&params=30_14_40.7_S_70_44_57.9_W_region:CL-CO_type:landmark
+    X05 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "X05")
+    lon, lat = X05.lon_lat()
+    np.testing.assert_allclose(lon[0], -70.749417, atol=1e-4, rtol=0)
+    np.testing.assert_allclose(lat[0], -30.244639, atol=1e-4, rtol=0)
+
+    # From: https://geohack.toolforge.org/geohack.php?pagename=V%C3%ADctor_M._Blanco_Telescope&params=30.16967_S_70.80653_W_
+    W84 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "W84")
+    lon, lat = W84.lon_lat()
+    np.testing.assert_allclose(lon[0], -70.80653, atol=1e-3, rtol=0)
+    np.testing.assert_allclose(lat[0], -30.169679, atol=1e-3, rtol=0)
+
+    # From: https://geohack.toolforge.org/geohack.php?params=32_22_48.6_S_20_48_38.1_E
+    M22 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "M22")
+    lon, lat = M22.lon_lat()
+    np.testing.assert_allclose(lon[0], 20.810583, atol=1e-4, rtol=0)
+    np.testing.assert_allclose(lat[0], -32.380167, atol=1e-4, rtol=0)
+
+    # From: https://geohack.toolforge.org/geohack.php?params=20_42_26.04_N_156_15_21.28_W
+    F51 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "F51")
+    lon, lat = F51.lon_lat()
+    np.testing.assert_allclose(lon[0], -156.255911, atol=1e-4, rtol=0)
+    np.testing.assert_allclose(lat[0], 20.707233, atol=1e-4, rtol=0)
+
+
+def test_ObservatoryParallaxCoeffiecients_timezone() -> None:
+    I41 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "I41")
+    assert I41.timezone() == "America/Los_Angeles"
+
+    X05 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "X05")
+    assert X05.timezone() == "America/Santiago"
+
+    W84 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "W84")
+    assert W84.timezone() == "America/Santiago"
+
+    M22 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "M22")
+    assert M22.timezone() == "Africa/Johannesburg"
+
+    F51 = OBSERVATORY_PARALLAX_COEFFICIENTS.select("code", "F51")
+    assert F51.timezone() == "Pacific/Honolulu"

src/adam_core/observers/tests/test_utils.py

@@ -0,0 +1,104 @@
+import pyarrow as pa
+import pyarrow.compute as pc
+
+from ...time import Timestamp
+from ..utils import calculate_observing_night
+
+
+def test_calculate_observing_night() -> None:
+    # Rubin Observatory is UTC-7
+    # Reasonable observating times would be +- 12 hours from local midnight
+    # or 7:00 to 17:00 UTC
+
+    # ZTF is UTC-8
+    # Reasonable observating times would be +- 12 hours from local midnight
+    # or 8:00 to 16:00 UTC
+
+    # M22 is UTC+2
+    # Reasonable observating times would be +- 12 hours from local midnight
+    # or 22:00 to 10:00 UTC
+
+    # 000 is UTC
+    # Reasonable observating times would be +- 12 hours from local midnight
+    # or 0:00 to 12:00 UTC
+
+    codes = pa.array(
+        [
+            "X05",
+            "X05",
+            "X05",
+            "X05",
+            "X05",
+            "I41",
+            "I41",
+            "I41",
+            "I41",
+            "I41",
+            "M22",
+            "M22",
+            "M22",
+            "M22",
+            "M22",
+            "000",
+            "000",
+            "000",
+            "000",
+            "000",
+        ]
+    )
+    times_utc = Timestamp.from_mjd(
+        [
+            # Rubin Observatory
+            59000 + 7 / 24 - 12 / 24,
+            59000 + 7 / 24 - 6 / 24,
+            59000 + 7 / 24,
+            59000 + 7 / 24 + 6 / 24,
+            59000 + 7 / 24 + 12 / 24,
+            # ZTF
+            59000 + 8 / 24 - 12 / 24,
+            59000 + 8 / 24 - 4 / 24,
+            59000 + 8 / 24,
+            59000 + 8 / 24 + 4 / 24,
+            59000 + 8 / 24 + 12 / 24,
+            # M22
+            59000 - 2 / 24 - 12 / 24,
+            59000 - 2 / 24 - 6 / 24,
+            59000 - 2 / 24,
+            59000 - 2 / 24 + 6 / 24,
+            59000 - 2 / 24 + 12 / 24,
+            # 000
+            59000 - 12 / 24,
+            59000 - 6 / 24,
+            59000,
+            59000 + 6 / 24,
+            59000 + 12 / 24,
+        ],
+        scale="utc",
+    )
+
+    observing_night = calculate_observing_night(codes, times_utc)
+    desired = pa.array(
+        [
+            58999,
+            58999,
+            58999,
+            58999,
+            59000,
+            58999,
+            58999,
+            58999,
+            58999,
+            59000,
+            58999,
+            58999,
+            58999,
+            58999,
+            59000,
+            58999,
+            58999,
+            58999,
+            58999,
+            59000,
+        ]
+    )
+    assert pc.all(pc.equal(observing_night, desired)).as_py()