Rewrite night day ratio algorithm (#313)

ecoscope/analysis/astronomy.py

@@ -1,13 +1,15 @@
 import warnings
+from datetime import datetime, timedelta
 
 import numpy as np
 import pandas as pd
 import pyproj
+import pytz
 
 try:
     import astroplan
-    import astropy.coordinates
-    import astropy.time
+    from astropy.coordinates import EarthLocation
+    from astropy.time import Time
 except ModuleNotFoundError:
     raise ModuleNotFoundError(
         'Missing optional dependencies required by this module. \
@@ -33,7 +35,7 @@ def to_EarthLocation(geometry):
         proj_from="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
         proj_to="+proj=geocent +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
     )
-    return astropy.coordinates.EarthLocation.from_geocentric(
+    return EarthLocation.from_geocentric(
         *trans.transform(xx=geometry.x, yy=geometry.y, zz=np.zeros(geometry.shape[0])), unit="m"
     )
 
@@ -44,70 +46,53 @@ def is_night(geometry, time):
         return astroplan.Observer(to_EarthLocation(geometry.centroid)).is_night(time)
 
 
-def get_daynight_ratio(traj, n_grid_points=150) -> pd.Series:
-    """
-    Parameters
-    ----------
-    n_grid_points : int (optional)
-        The number of grid points on which to search for the horizon
-        crossings of the target over a 24-hour period, default is 150 which
-        yields rise time precisions better than one minute.
-        https://github.com/astropy/astroplan/pull/424
-    Returns
-    -------
-    pd.Series:
-        Daynight ratio for each unique individual subject in the grouby_col column.
-    """
+def sun_time(date, geometry):
+    midnight = Time(
+        datetime(date.year, date.month, date.day) + timedelta(seconds=1), scale="utc"
+    )  # add 1 second shift to avoid leap_second_strict warning
+    observer = astroplan.Observer(location=EarthLocation(lon=geometry.centroid.x, lat=geometry.centroid.y))
+    sunrise = observer.sun_rise_time(midnight, which="next", n_grid_points=150).to_datetime(timezone=pytz.UTC)
+    sunset = observer.sun_set_time(midnight, which="next", n_grid_points=150).to_datetime(timezone=pytz.UTC)
+    return pd.Series({"sunrise": sunrise, "sunset": sunset})
+
+
+def calculate_day_night_distance(date, segment_start, segment_end, dist_meters, daily_summary):
+    sunrise = daily_summary.loc[date, "sunrise"]
+    sunset = daily_summary.loc[date, "sunset"]
+
+    if segment_start < sunset and segment_end > sunset:  # start in day and end in night
+        day_percent = (sunset - segment_start) / (segment_end - segment_start)
+    elif segment_start < sunrise and segment_end > sunrise:  # start in night and end in day
+        day_percent = (segment_end - sunrise) / (segment_end - segment_start)
+    elif sunrise < sunset:
+        if segment_end < sunrise or segment_start > sunset:  # all night
+            day_percent = 0
+        elif segment_start >= sunrise and segment_end <= sunset:  # all day
+            day_percent = 1
+    else:  # sunrise >= sunset
+        if segment_end < sunset or segment_start > sunrise:  # all day
+            day_percent = 1
+        elif segment_start >= sunset and segment_end <= sunrise:  # all night
+            day_percent = 0
+
+    daily_summary.loc[date, "day_distance"] += day_percent * dist_meters
+    daily_summary.loc[date, "night_distance"] += (1 - day_percent) * dist_meters
+
+
+def get_nightday_ratio(gdf):
+    gdf["date"] = pd.to_datetime(gdf["segment_start"]).dt.date
+
+    daily_summary = gdf.groupby("date").first()["geometry"].reset_index()
+    daily_summary[["sunrise", "sunset"]] = daily_summary.apply(lambda x: sun_time(x.date, x.geometry), axis=1)
+    daily_summary["day_distance"] = 0.0
+    daily_summary["night_distance"] = 0.0
+    daily_summary = daily_summary.set_index("date")
+
+    gdf.apply(
+        lambda x: calculate_day_night_distance(x.date, x.segment_start, x.segment_end, x.dist_meters, daily_summary),
+        axis=1,
+    )
 
-    locations = to_EarthLocation(traj.geometry.to_crs(crs=traj.estimate_utm_crs()).centroid)
-
-    observer = astroplan.Observer(location=locations)
-    is_night_start = observer.is_night(traj.segment_start)
-    is_night_end = observer.is_night(traj.segment_end)
-
-    # Night -> Night
-    night_distance = traj.dist_meters.loc[is_night_start & is_night_end].sum()
-
-    # Day -> Day
-    day_distance = traj.dist_meters.loc[~is_night_start & ~is_night_end].sum()
-
-    # Night -> Day
-    night_day_mask = is_night_start & ~is_night_end
-    night_day_df = traj.loc[night_day_mask, ["segment_start", "dist_meters", "timespan_seconds"]]
-    i = (
-        pd.to_datetime(
-            astroplan.Observer(location=locations[night_day_mask])
-            .sun_rise_time(
-                astropy.time.Time(night_day_df.segment_start),
-                n_grid_points=n_grid_points,
-                which="next",
-            )
-            .datetime,
-            utc=True,
-        )
-        - night_day_df.segment_start
-    ).dt.total_seconds() / night_day_df.timespan_seconds
-
-    night_distance += (night_day_df.dist_meters * i).sum()
-    day_distance += ((1 - i) * night_day_df.dist_meters).sum()
-
-    # Day -> Night
-    day_night_mask = ~is_night_start & is_night_end
-    day_night_df = traj.loc[day_night_mask, ["segment_start", "dist_meters", "timespan_seconds"]]
-    i = (
-        pd.to_datetime(
-            astroplan.Observer(location=locations[day_night_mask])
-            .sun_set_time(
-                astropy.time.Time(day_night_df.segment_start),
-                n_grid_points=n_grid_points,
-                which="next",
-            )
-            .datetime,
-            utc=True,
-        )
-        - day_night_df.segment_start
-    ).dt.total_seconds() / day_night_df.timespan_seconds
-    day_distance += (day_night_df.dist_meters * i).sum()
-    night_distance += ((1 - i) * day_night_df.dist_meters).sum()
-
-    return night_distance / day_distance
+    daily_summary["night_day_ratio"] = daily_summary["night_distance"] / daily_summary["day_distance"]
+    mean_night_day_ratio = daily_summary["night_day_ratio"].replace([np.inf, -np.inf], np.nan).dropna().mean()
+    return mean_night_day_ratio

tests/test_astronomy.py

@@ -1,5 +1,8 @@
+from datetime import datetime
+
 import pandas as pd
 import pyproj
+import pytest
 
 from ecoscope.analysis import astronomy
 from ecoscope.base import Trajectory
@@ -33,16 +36,121 @@ def test_is_night(movebank_relocations):
     assert subset["is_night"].values.tolist() == [True, True, False]
 
 
-def test_daynight_ratio(movebank_relocations):
+def test_nightday_ratio(movebank_relocations):
     trajectory = Trajectory.from_relocations(movebank_relocations)
     expected = pd.Series(
-        [
-            0.451912,
-            1.523379,
-        ],
+        [0.45905845612291696, 2.0019632541788472],
         index=pd.Index(["Habiba", "Salif Keita"], name="groupby_col"),
     )
     pd.testing.assert_series_equal(
-        trajectory.groupby("groupby_col")[trajectory.columns].apply(astronomy.get_daynight_ratio, include_groups=False),
+        trajectory.groupby("groupby_col")[trajectory.columns].apply(astronomy.get_nightday_ratio, include_groups=False),
         expected,
     )
+
+
+@pytest.fixture
+def daily_summary():
+    """Create a sample daily summary DataFrame for testing."""
+    date = datetime(2024, 1, 1)
+    df = pd.DataFrame(
+        {
+            "sunrise": [datetime(2024, 1, 1, 6, 0)],
+            "sunset": [datetime(2024, 1, 1, 18, 0)],
+            "day_distance": [0.0],
+            "night_distance": [0.0],
+        },
+        index=[date],
+    )
+    return df
+
+
+def test_all_night_before_sunrise(daily_summary):
+    """Test segment entirely before sunrise."""
+    date = datetime(2024, 1, 1)
+    astronomy.calculate_day_night_distance(
+        date, datetime(2024, 1, 1, 2, 0), datetime(2024, 1, 1, 4, 0), 1000, daily_summary
+    )  # 2:00 AM  # 4:00 AM
+    assert daily_summary.loc[date, "night_distance"] == 1000
+    assert daily_summary.loc[date, "day_distance"] == 0
+
+
+def test_all_night_after_sunset(daily_summary):
+    """Test segment entirely after sunset."""
+    date = datetime(2024, 1, 1)
+    astronomy.calculate_day_night_distance(
+        date, datetime(2024, 1, 1, 20, 0), datetime(2024, 1, 1, 22, 0), 1000, daily_summary
+    )  # 8:00 PM  # 10:00 PM
+    assert daily_summary.loc[date, "night_distance"] == 1000
+    assert daily_summary.loc[date, "day_distance"] == 0
+
+
+def test_all_day(daily_summary):
+    """Test segment entirely during daylight hours."""
+    date = datetime(2024, 1, 1)
+    astronomy.calculate_day_night_distance(
+        date, datetime(2024, 1, 1, 10, 0), datetime(2024, 1, 1, 14, 0), 1000, daily_summary
+    )  # 10:00 AM  # 2:00 PM
+    assert daily_summary.loc[date, "day_distance"] == 1000
+    assert daily_summary.loc[date, "night_distance"] == 0
+
+
+def test_day_to_night_transition(daily_summary):
+    """Test segment starting in day and ending in night."""
+    date = datetime(2024, 1, 1)
+    astronomy.calculate_day_night_distance(
+        date, datetime(2024, 1, 1, 17, 0), datetime(2024, 1, 1, 19, 0), 1000, daily_summary
+    )  # 5:00 PM  # 7:00 PM
+    # Segment spans 2 hours, with 1 hour in day and 1 hour in night
+    assert daily_summary.loc[date, "day_distance"] == 500
+    assert daily_summary.loc[date, "night_distance"] == 500
+
+
+def test_night_to_day_transition(daily_summary):
+    """Test segment starting in night and ending in day."""
+    date = datetime(2024, 1, 1)
+    astronomy.calculate_day_night_distance(
+        date, datetime(2024, 1, 1, 5, 0), datetime(2024, 1, 1, 7, 0), 1000, daily_summary
+    )  # 5:00 AM  # 7:00 AM
+    # Segment spans 2 hours, with 1 hour in night and 1 hour in day
+    assert daily_summary.loc[date, "day_distance"] == 500
+    assert daily_summary.loc[date, "night_distance"] == 500
+
+
+@pytest.fixture
+def inverted_daily_summary():
+    """Create a daily summary where sunrise is after sunset in UTC."""
+    date = datetime(2024, 1, 1)
+    df = pd.DataFrame(
+        {
+            "sunrise": [datetime(2024, 1, 1, 18, 0)],  # 6:00 PM
+            "sunset": [datetime(2024, 1, 1, 6, 0)],  # 6:00 AM
+            "day_distance": [0.0],
+            "night_distance": [0.0],
+        },
+        index=[date],
+    )
+    return df
+
+
+def test_inverted_all_day(inverted_daily_summary):
+    """Test all-day segment when sunrise is after sunset."""
+    date = datetime(2024, 1, 1)
+    astronomy.calculate_day_night_distance(
+        date, datetime(2024, 1, 1, 4, 0), datetime(2024, 1, 1, 5, 0), 1000, inverted_daily_summary  # 4:00 AM  # 5:00 AM
+    )
+    assert inverted_daily_summary.loc[date, "day_distance"] == 1000
+    assert inverted_daily_summary.loc[date, "night_distance"] == 0
+
+
+def test_inverted_all_night(inverted_daily_summary):
+    """Test all-night segment when sunrise is after sunset."""
+    date = datetime(2024, 1, 1)
+    astronomy.calculate_day_night_distance(
+        date,
+        datetime(2024, 1, 1, 7, 0),
+        datetime(2024, 1, 1, 17, 0),
+        1000,
+        inverted_daily_summary,  # 7:00 AM  # 5:00 PM
+    )
+    assert inverted_daily_summary.loc[date, "night_distance"] == 1000
+    assert inverted_daily_summary.loc[date, "day_distance"] == 0