Support decimal-year epochs in GeodePy transformations

geodepy/constants.py

@@ -5,13 +5,35 @@
 Constants Module
 """
 
-from math import sqrt
-from datetime import date
+from datetime import date, datetime, time
+from math import floor, isfinite, sqrt
 
 c_vac = 299792.458
 k_0 = 0.9996
 
 
+# Epoch conversion
+def _epoch_datetime(value):
+    """Convert a calendar date or decimal epoch to a datetime.
+
+    Decimal fractions count through the actual Gregorian calendar year:
+    2024.5 is 2024-07-02 00:00; 2025.5 is 2025-07-02 12:00.
+    The transformation's existing 365.25-day rate denominator is unchanged.
+    """
+    if type(value) == datetime:
+        return value
+    if type(value) == date:
+        return datetime.combine(value, time.min)
+    if type(value) not in (int, float) or not isfinite(value):
+        raise ValueError("epoch must be a date or finite decimal year")
+    year = floor(value)
+    if not 1 <= year <= 9998:
+        raise ValueError("decimal year outside supported calendar range")
+    start = datetime(year, 1, 1)
+    duration = datetime(year + 1, 1, 1) - start
+    return start + duration * (value - year)
+
+
 # Ellipsoid constants
 class Ellipsoid(object):
     def __init__(self, semimaj, inversef):
@@ -237,59 +259,61 @@ def __add__(self, other):
         Change the transformation epoch to a specified date. Advances all
         transformation parameters by their respective rates of change.
 
-        :param other: datetime.date object
+        :param other: datetime.date object or decimal epoch (YYYY.fraction)
         :return: transformation object with parameters and ref epoch moved to
                  the specified date
         """
-        if type(other) == date:
+        if type(other) is date and type(self.ref_epoch) is date:
+            # Preserve existing calendar-date behaviour exactly.
             timediff = (other - self.ref_epoch).days / 365.25
-
-            if type(self.tf_sd) == TransformationSD:
-                self.tf_sd.sd_tx = (
-                    self.tf_sd.sd_tx**2 + (self.tf_sd.sd_d_tx * timediff) ** 2
-                ) ** 0.5
-                self.tf_sd.sd_ty = (
-                    self.tf_sd.sd_ty**2 + (self.tf_sd.sd_d_ty * timediff) ** 2
-                ) ** 0.5
-                self.tf_sd.sd_tz = (
-                    self.tf_sd.sd_tz**2 + (self.tf_sd.sd_d_tz * timediff) ** 2
-                ) ** 0.5
-                self.tf_sd.sd_sc = (
-                    self.tf_sd.sd_sc**2 + (self.tf_sd.sd_d_sc * timediff) ** 2
-                ) ** 0.5
-                self.tf_sd.sd_rx = (
-                    self.tf_sd.sd_rx**2 + (self.tf_sd.sd_d_rx * timediff) ** 2
-                ) ** 0.5
-                self.tf_sd.sd_ry = (
-                    self.tf_sd.sd_ry**2 + (self.tf_sd.sd_d_ry * timediff) ** 2
-                ) ** 0.5
-                self.tf_sd.sd_rz = (
-                    self.tf_sd.sd_rz**2 + (self.tf_sd.sd_d_rz * timediff) ** 2
-                ) ** 0.5
-
-            return Transformation(
-                self.to_datum,
-                self.from_datum,
-                other,
-                round(self.tx + (self.d_tx * timediff), 8),
-                round(self.ty + (self.d_ty * timediff), 8),
-                round(self.tz + (self.d_tz * timediff), 8),
-                round(self.sc + (self.d_sc * timediff), 8),
-                round(self.rx + (self.d_rx * timediff), 8),
-                round(self.ry + (self.d_ry * timediff), 8),
-                round(self.rz + (self.d_rz * timediff), 8),
-                self.d_tx,
-                self.d_ty,
-                self.d_tz,
-                self.d_sc,
-                self.d_rx,
-                self.d_ry,
-                self.d_rz,
-                self.tf_sd,
-            )
         else:
-            ValueError("supports adding datetime.date objects only")
+            timediff = (
+                _epoch_datetime(other) - _epoch_datetime(self.ref_epoch)
+            ).total_seconds() / (86400 * 365.25)
+
+        if type(self.tf_sd) == TransformationSD:
+            self.tf_sd.sd_tx = (
+                self.tf_sd.sd_tx**2 + (self.tf_sd.sd_d_tx * timediff) ** 2
+            ) ** 0.5
+            self.tf_sd.sd_ty = (
+                self.tf_sd.sd_ty**2 + (self.tf_sd.sd_d_ty * timediff) ** 2
+            ) ** 0.5
+            self.tf_sd.sd_tz = (
+                self.tf_sd.sd_tz**2 + (self.tf_sd.sd_d_tz * timediff) ** 2
+            ) ** 0.5
+            self.tf_sd.sd_sc = (
+                self.tf_sd.sd_sc**2 + (self.tf_sd.sd_d_sc * timediff) ** 2
+            ) ** 0.5
+            self.tf_sd.sd_rx = (
+                self.tf_sd.sd_rx**2 + (self.tf_sd.sd_d_rx * timediff) ** 2
+            ) ** 0.5
+            self.tf_sd.sd_ry = (
+                self.tf_sd.sd_ry**2 + (self.tf_sd.sd_d_ry * timediff) ** 2
+            ) ** 0.5
+            self.tf_sd.sd_rz = (
+                self.tf_sd.sd_rz**2 + (self.tf_sd.sd_d_rz * timediff) ** 2
+            ) ** 0.5
 
+        return Transformation(
+            self.to_datum,
+            self.from_datum,
+            other,
+            round(self.tx + (self.d_tx * timediff), 8),
+            round(self.ty + (self.d_ty * timediff), 8),
+            round(self.tz + (self.d_tz * timediff), 8),
+            round(self.sc + (self.d_sc * timediff), 8),
+            round(self.rx + (self.d_rx * timediff), 8),
+            round(self.ry + (self.d_ry * timediff), 8),
+            round(self.rz + (self.d_rz * timediff), 8),
+            self.d_tx,
+            self.d_ty,
+            self.d_tz,
+            self.d_sc,
+            self.d_rx,
+            self.d_ry,
+            self.d_rz,
+            self.tf_sd,
+        )
 
 class TransformationSD(object):
     def __init__(

geodepy/tests/test_decimal_epoch.py

@@ -0,0 +1,132 @@
+"""Regression tests for decimal-year epoch support (GeodePy #104)."""
+
+from datetime import date
+
+import pytest
+
+from geodepy.constants import _epoch_datetime, itrf2014_to_gda2020
+from geodepy.transform import (
+    conform14,
+    plate_motion_transformation,
+    transform_atrf2014_to_gda2020,
+)
+
+COORD = (-4052052.7373, 4212835.9835, -2545104.5867)
+MODEL = itrf2014_to_gda2020
+
+
+@pytest.mark.parametrize("year", [1994, 2000, 2018, 2020, 2024, 2025, 2026])
+def test_integer_decimal_year_matches_january_first(year):
+    assert conform14(*COORD, float(year), MODEL) == conform14(
+        *COORD, date(year, 1, 1), MODEL
+    )
+
+
+@pytest.mark.parametrize(
+    "from_year,to_year", [(2018, 2025), (2025, 2018), (2020, 2020), (2024, 2026)]
+)
+def test_plate_motion_integer_years_match_dates(from_year, to_year):
+    expected = plate_motion_transformation(
+        *COORD, date(from_year, 1, 1), date(to_year, 1, 1), MODEL
+    )
+    assert plate_motion_transformation(
+        *COORD, float(from_year), float(to_year), MODEL
+    ) == expected
+
+
+@pytest.mark.parametrize(
+    "from_epoch,to_epoch",
+    [(2018.0, date(2025, 1, 1)), (date(2018, 1, 1), 2025.0)],
+)
+def test_mixed_input_types(from_epoch, to_epoch):
+    expected = plate_motion_transformation(
+        *COORD, date(2018, 1, 1), date(2025, 1, 1), MODEL
+    )
+    assert plate_motion_transformation(*COORD, from_epoch, to_epoch, MODEL) == expected
+
+
+@pytest.mark.parametrize(
+    "year,fraction,expected",
+    [
+        (2024, 0.0, (2024, 1, 1, 0)),
+        (2024, 0.5, (2024, 7, 2, 0)),
+        (2025, 0.5, (2025, 7, 2, 12)),
+        (2026, 0.0, (2026, 1, 1, 0)),
+        (2026, 0.5, (2026, 7, 2, 12)),
+    ],
+)
+def test_decimal_calendar_fraction(year, fraction, expected):
+    observed = _epoch_datetime(year + fraction)
+    assert (observed.year, observed.month, observed.day, observed.hour) == expected
+
+
+def test_leap_halfway_matches_same_calendar_date():
+    assert conform14(*COORD, 2024.5, MODEL) == conform14(
+        *COORD, date(2024, 7, 2), MODEL
+    )
+
+
+def test_fractional_year_affects_transform():
+    jan = conform14(*COORD, 2025.0, MODEL)
+    mid = conform14(*COORD, 2025.5, MODEL)
+    next_year = conform14(*COORD, 2026.0, MODEL)
+    for a, b, c in zip(jan[:3], mid[:3], next_year[:3]):
+        assert abs(a - c) > 0
+        assert min(a, c) <= b <= max(a, c)
+
+
+def test_fractional_plate_motion_reverse_round_trip():
+    endpoint = plate_motion_transformation(*COORD, 2018.25, 2025.75, MODEL)
+    reverse = plate_motion_transformation(*endpoint[:3], 2025.75, 2018.25, MODEL)
+    for actual, original in zip(reverse[:3], COORD):
+        assert abs(actual - original) < 0.002
+
+
+def test_wrapper_supports_decimal_year():
+    assert transform_atrf2014_to_gda2020(*COORD, 2018.0) == (
+        transform_atrf2014_to_gda2020(*COORD, date(2018, 1, 1))
+    )
+
+
+@pytest.mark.parametrize(
+    "bad", [float("nan"), float("inf"), float("-inf"), 0.5, -2026.1, 10000.0]
+)
+def test_invalid_year_fails(bad):
+    with pytest.raises(ValueError):
+        conform14(*COORD, bad, MODEL)
+
+
+@pytest.mark.parametrize("bad", [True, False, "2025.5", None])
+def test_unsupported_year_type_fails(bad):
+    with pytest.raises(ValueError):
+        conform14(*COORD, bad, MODEL)
+
+
+@pytest.mark.parametrize(
+    "invalid",
+    [float("nan"), float("inf"), float("-inf"), True, "2025.5", None, 0.1, 10000.0],
+)
+@pytest.mark.parametrize("slot", [0, 1])
+def test_plate_motion_invalid_decimal_years(invalid, slot):
+    epochs = [2018.25, 2025.75]
+    epochs[slot] = invalid
+    with pytest.raises(ValueError):
+        plate_motion_transformation(*COORD, *epochs, MODEL)
+
+
+def test_invalid_transform_addition_fails_loudly():
+    with pytest.raises(ValueError):
+        MODEL + float("nan")
+    with pytest.raises(ValueError):
+        MODEL + "2018.0"
+
+
+def test_ref_epoch_preserved_on_date_path():
+    result = MODEL + date(2026, 5, 11)
+    assert result.ref_epoch == date(2026, 5, 11)
+    assert type(result.ref_epoch) is date
+
+
+def test_ref_epoch_is_numeric_on_numeric_path():
+    result = MODEL + 2026.25
+    assert result.ref_epoch == 2026.25

geodepy/transform.py

@@ -17,6 +17,7 @@
 from geodepy.constants import (
     Transformation,
     TransformationSD,
+    _epoch_datetime,
     atrf2014_to_gda2020,
     gda94_to_gda2020,
 )
@@ -130,15 +131,15 @@ def conform14(x, y, z, to_epoch, trans, vcv=None):
     :param x: Cartesian X (m)
     :param y: Cartesian Y (m)
     :param z: Cartesian Z (m)
-    :param to_epoch: Epoch co-ordinate transformation is performed at (datetime.date Object)
+    :param to_epoch: Target epoch (datetime.date or decimal year, e.g. 2025.5)
     :param trans: Transformation Object
     :param vcv: Optional 3*3 numpy array in Cartesian units to propagate tf uncertainty
     :return: Cartesian X, Y, Z co-ordinates and vcv matrix transformed using Transformation parameters at desired epoch
     """
     if type(trans) != Transformation:
         raise ValueError("trans must be a Transformation Object")
-    if type(to_epoch) != datetime.date:
-        raise ValueError("to_epoch must be a datetime.date Object")
+    if type(to_epoch) not in (datetime.date, int, float):
+        raise ValueError("to_epoch must be a datetime.date or decimal year")
     # Calculate 7 Parameters from 14 Parameter Transformation Object
     timetrans = trans + to_epoch
 
@@ -153,24 +154,26 @@ def plate_motion_transformation(x, y, z, from_epoch, to_epoch, plate_motion, vcv
     :param x: Cartesian X (m)
     :param y: Cartesian Y (m)
     :param z: Cartesian Z (m)
-    :param from_epoch: Epoch the co-ordinate transformation is from (datetime.date Object)
-    :param to_epoch: Epoch the co-ordinate transformation is to (datetime.date Object)
+    :param from_epoch: Source epoch (datetime.date or decimal year)
+    :param to_epoch: Target epoch (datetime.date or decimal year)
     :param plate_motion: Plate motion model for transformation
     :param vcv: Optional 3*3 numpy array in Cartesian units to propagate tf uncertainty
     :return: Cartesian X, Y, Z co-ordinates and vcv matrix transformed using plate motion to desired epoch
     """
-    if type(to_epoch) != datetime.date:
-        raise ValueError("to_epoch must be a datetime.date Object")
-    if type(from_epoch) != datetime.date:
-        raise ValueError("from_epoch must be a datetime.date Object")
+    if type(to_epoch) not in (datetime.date, int, float):
+        raise ValueError("to_epoch must be a datetime.date or decimal year")
+    if type(from_epoch) not in (datetime.date, int, float):
+        raise ValueError("from_epoch must be a datetime.date or decimal year")
     if type(plate_motion) != Transformation:
         raise ValueError("plate_motion must be a Transformation Object")
 
-    #calculate number of years to be moved
-    timediff= to_epoch - from_epoch
-
-    #calculate epoch needed for plate motion
-    change_epoch = plate_motion.ref_epoch - timediff
+    # Retain exact existing day arithmetic for calendar-date inputs.
+    if type(from_epoch) is datetime.date and type(to_epoch) is datetime.date:
+        timediff = to_epoch - from_epoch
+        change_epoch = plate_motion.ref_epoch - timediff
+    else:
+        timediff = _epoch_datetime(to_epoch) - _epoch_datetime(from_epoch)
+        change_epoch = _epoch_datetime(plate_motion.ref_epoch) - timediff
 
     # Calculate 7 Parameters from 14 Parameter Transformation Object
     timetrans = plate_motion + change_epoch