From 201af4c6ef8393519e7efa57666c042649792ff3 Mon Sep 17 00:00:00 2001 From: Yun-Wu Date: Fri, 15 Nov 2024 20:18:35 +0800 Subject: [PATCH] Rewrite night day ratio algorithm (#313) --- ecoscope/analysis/astronomy.py | 123 +++++++++++++++------------------ tests/test_astronomy.py | 120 ++++++++++++++++++++++++++++++-- 2 files changed, 168 insertions(+), 75 deletions(-) diff --git a/ecoscope/analysis/astronomy.py b/ecoscope/analysis/astronomy.py index 825fb23d..7b7aafb4 100644 --- a/ecoscope/analysis/astronomy.py +++ b/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 diff --git a/tests/test_astronomy.py b/tests/test_astronomy.py index 50f53e75..66f881bd 100644 --- a/tests/test_astronomy.py +++ b/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