diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index 20c7ca1..30c4e71 100644 --- a/pyorbital/orbital.py +++ b/pyorbital/orbital.py @@ -68,60 +68,77 @@ class OrbitalError(Exception): pass -def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): - """Calculate observers look angle to a satellite. - - http://celestrak.com/columns/v02n02/ - - :utc_time: Observation time (datetime object) - :lon: Longitude of observer position on ground in degrees east - :lat: Latitude of observer position on ground in degrees north - :alt: Altitude above sea-level (geoid) of observer position on ground in km - - :return: (Azimuth, Elevation) - """ - (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = astronomy.observer_position( - utc_time, sat_lon, sat_lat, sat_alt) - - (opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \ - astronomy.observer_position(utc_time, lon, lat, alt) - - lon = np.deg2rad(lon) - lat = np.deg2rad(lat) - - theta = (astronomy.gmst(utc_time) + lon) % (2 * np.pi) - - rx = pos_x - opos_x - ry = pos_y - opos_y - rz = pos_z - opos_z - +def ecef_to_topocentric(rx, ry, rz, lat, lon, theta): + """Convert ECEF vector to topocentric-horizon coordinates.""" sin_lat = np.sin(lat) cos_lat = np.cos(lat) sin_theta = np.sin(theta) cos_theta = np.cos(theta) - top_s = sin_lat * cos_theta * rx + \ - sin_lat * sin_theta * ry - cos_lat * rz + # Transform ECEF coordinates to topocentric-horizon coordinates + top_s = sin_lat * cos_theta * rx + sin_lat * sin_theta * ry - cos_lat * rz top_e = -sin_theta * rx + cos_theta * ry - top_z = cos_lat * cos_theta * rx + \ - cos_lat * sin_theta * ry + sin_lat * rz + top_z = cos_lat * cos_theta * rx + cos_lat * sin_theta * ry + sin_lat * rz + + return top_s, top_e, top_z +def compute_azimuth_elevation(top_s, top_e, top_z, rg_): + """Compute azimuth and elevation from topocentric coordinates.""" # Azimuth is undefined when elevation is 90 degrees, 180 (pi) will be returned. az_ = np.arctan2(-top_e, top_s) + np.pi az_ = np.mod(az_, 2 * np.pi) # Needed on some platforms - rg_ = np.sqrt(rx * rx + ry * ry + rz * rz) - - top_z_divided_by_rg_ = top_z / rg_ - # Due to rounding top_z can be larger than rg_ (when el_ ~ 90). - top_z_divided_by_rg_ = top_z_divided_by_rg_.clip(max=1) + top_z_divided_by_rg_ = np.clip(top_z / rg_, -1, 1) el_ = np.arcsin(top_z_divided_by_rg_) return np.rad2deg(az_), np.rad2deg(el_) +def _look_from_ecef(pos, opos, lon_rad, lat_rad, theta): + """Compute azimuth/elevation from two ECEF positions.""" + rx = pos[0] - opos[0] + ry = pos[1] - opos[1] + rz = pos[2] - opos[2] + rg_ = np.sqrt(rx * rx + ry * ry + rz * rz) -class Orbital(object): + top_s, top_e, top_z = ecef_to_topocentric(rx, ry, rz, lat_rad, lon_rad, theta) + return compute_azimuth_elevation(top_s, top_e, top_z, rg_) + +def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): + """Calculate observer's look angle to a satellite. + + http://celestrak.com/columns/v02n02/ + + :param sat_lon: Satellite longitude in degrees east + :param sat_lat: Satellite latitude in degrees north + :param sat_alt: Satellite altitude in km + :param utc_time: Observation time (datetime object) + :param lon: Observer longitude in degrees east + :param lat: Observer latitude in degrees north + :param alt: Observer altitude in km + :return: (Azimuth, Elevation) in degrees + """ + (pos_x, pos_y, pos_z), _ = astronomy.observer_position( + utc_time, sat_lon, sat_lat, sat_alt + ) + (opos_x, opos_y, opos_z), _ = astronomy.observer_position( + utc_time, lon, lat, alt + ) + + lon_rad = np.deg2rad(lon) + lat_rad = np.deg2rad(lat) + theta = (astronomy.gmst(utc_time) + lon_rad) % (2 * np.pi) + + return _look_from_ecef( + (pos_x, pos_y, pos_z), + (opos_x, opos_y, opos_z), + lon_rad, + lat_rad, + theta + ) + + +class Orbital: """Class for orbital computations. The *satellite* parameter is the name of the satellite to work on and is @@ -144,39 +161,60 @@ def __str__(self): """Print the Orbital object state.""" return self.satellite_name + " " + str(self.tle) + def _find_last_node_time(self, utc_time, is_ascending=True): + """Find the last node crossing time (ascending or descending) before utc_time.""" + time_ref = np.datetime64(_get_tz_unaware_utctime(utc_time)) + mean_motion = self.tle.mean_motion + orbit_period_min = XMNPDA / mean_motion + + def node_func(minutes: float) -> float: + """Returns the satellite Z-position (in km) at time_ref + minutes.""" + time_f = time_ref + np.timedelta64(int(minutes * 60), "s") + pos, _ = self.get_position(time_f, normalize=False) + return pos[2] + + def find_bracket(func, start_min, step_min, max_back_min): + """Find a valid bracket where func crosses zero (sign change).""" + t_high = start_min + t_low = t_high - step_min + while abs(t_low) <= max_back_min: + f_high = func(t_high) + f_low = func(t_low) + if f_high * f_low < 0: + return t_low, t_high + t_high = t_low + t_low -= step_min + raise ValueError("Could not find suitable bracket for node crossing.") + + t_low_min, t_high_min = find_bracket( + node_func, + start_min=0.0, + step_min=1.0, + max_back_min=orbit_period_min * 2, + ) + + root_min = _get_root(node_func, t_low_min, t_high_min, tol=0.0001) + t_node = time_ref + np.timedelta64(int(root_min * 60), "s") + _, vel = self.get_position(t_node, normalize=False) + + if is_ascending and vel[2] < 0: + t_node -= np.timedelta64(int(orbit_period_min / 2.0 * 60), "s") + elif not is_ascending and vel[2] > 0: + t_node -= np.timedelta64(int(orbit_period_min / 2.0 * 60), "s") + + return t_node + def get_last_an_time(self, utc_time): """Calculate time of last ascending node relative to the specified time.""" - # Propagate backwards to ascending node - dt = np.timedelta64(10, "m") - - t_old = np.datetime64(_get_tz_unaware_utctime(utc_time)) - t_new = t_old - dt - pos0, vel0 = self.get_position(t_old, normalize=False) - pos1, vel1 = self.get_position(t_new, normalize=False) - while not (pos0[2] > 0 and pos1[2] < 0): - pos0 = pos1 - t_old = t_new - t_new = t_old - dt - pos1, vel1 = self.get_position(t_new, normalize=False) - - # Return if z within 1 km of an - if np.abs(pos0[2]) < 1: - return t_old - elif np.abs(pos1[2]) < 1: - return t_new - - # Bisect to z within 1 km - while np.abs(pos1[2]) > 1: - # pos0, vel0 = pos1, vel1 - dt = (t_old - t_new) / 2 - t_mid = t_old - dt - pos1, vel1 = self.get_position(t_mid, normalize=False) - if pos1[2] > 0: - t_old = t_mid - else: - t_new = t_mid + t_an = self._find_last_node_time(utc_time, is_ascending=True) + logger.debug(f"Ascending Node crossing time: {t_an}") + return t_an - return t_mid + def get_last_dn_time(self, utc_time): + """Calculate time of last descending node relative to the specified time.""" + t_dn = self._find_last_node_time(utc_time, is_ascending=False) + logger.debug(f"Descending Node crossing time: {t_dn}") + return t_dn def get_position(self, utc_time, normalize=True): """Get the cartesian position and velocity from the satellite.""" @@ -216,13 +254,70 @@ def get_lonlatalt(self, utc_time): alt *= A return np.rad2deg(lon), np.rad2deg(lat), alt - def find_aos(self, utc_time, lon, lat): - """Find AOS.""" - pass + def _find_single_crossing( + self, + utc_time: dt.datetime, + lon: float, + lat: float, + alt: float, + horizon: float, + is_aos: bool, + max_search_min: float = 1440.0 + ) -> dt.datetime | None: + """Internal helper to find the single next horizon crossing time (AOS or AOL).""" + elev_func = partial(self._elevation, utc_time, lon, lat, alt, horizon) - def find_aol(self, utc_time, lon, lat): - """Find AOL.""" - pass + t_min = 0.0 # Start time in minutes from utc_time + t_step = 5.0 # Initial search step + t_curr = t_min + + # Iterate forward until a sign change is found + while t_curr < max_search_min: + t_next = t_curr + t_step + + elev_curr = elev_func(t_curr) + elev_next = elev_func(t_next) + + # Check for a sign change (crossing the horizon) + if elev_curr * elev_next < 0: + t_low, t_high = sorted([t_curr, t_next]) + + # Check for rising (AOS) or setting (AOL) + # If elev_next > elev_curr, the satellite is rising (positive derivative) + is_rising = elev_next > elev_curr + + if (is_aos and is_rising) or (not is_aos and not is_rising): + root_min = optimize.brentq( + elev_func, t_low, t_high, xtol=0.00001 # High precision + ) + + return utc_time + dt.timedelta(minutes=root_min) + + t_curr = t_next + + return None # No crossing found within the search limit + + def find_aos( + self, + utc_time: dt.datetime, + lon: float, + lat: float, + alt: float = 0, + horizon: float = 0 + ) -> dt.datetime | None: + """Find the time of the next Acquisition of Signal (AOS) after utc_time (while rising).""" + return self._find_single_crossing(utc_time, lon, lat, alt, horizon, is_aos=True) + + def find_aol( + self, + utc_time: dt.datetime, + lon: float, + lat: float, + alt: float = 0, + horizon: float = 0 + ) -> dt.datetime | None: + """Find the time of the next Acquisition of Signal (AOL) after utc_time (while setting).""" + return self._find_single_crossing(utc_time, lon, lat, alt, horizon, is_aos=False) def get_observer_look(self, utc_time, lon, lat, alt): """Calculate observers look angle to a satellite. @@ -237,41 +332,22 @@ def get_observer_look(self, utc_time, lon, lat, alt): Return: (Azimuth, Elevation) """ - utc_time = dt2np(utc_time) - (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = self.get_position( - utc_time, normalize=False) - (opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \ - astronomy.observer_position(utc_time, lon, lat, alt) - - lon = np.deg2rad(lon) - lat = np.deg2rad(lat) + time_ref = dt2np(utc_time) - theta = (astronomy.gmst(utc_time) + lon) % (2 * np.pi) + (pos_x, pos_y, pos_z), _ = self.get_position(time_ref, normalize=False) + (opos_x, opos_y, opos_z), _ = astronomy.observer_position(time_ref, lon, lat, alt) - rx = pos_x - opos_x - ry = pos_y - opos_y - rz = pos_z - opos_z + lon_rad = np.deg2rad(lon) + lat_rad = np.deg2rad(lat) + theta = (astronomy.gmst(time_ref) + lon_rad) % (2 * np.pi) - sin_lat = np.sin(lat) - cos_lat = np.cos(lat) - sin_theta = np.sin(theta) - cos_theta = np.cos(theta) - - top_s = sin_lat * cos_theta * rx + \ - sin_lat * sin_theta * ry - cos_lat * rz - top_e = -sin_theta * rx + cos_theta * ry - top_z = cos_lat * cos_theta * rx + \ - cos_lat * sin_theta * ry + sin_lat * rz - - az_ = np.arctan(-top_e / top_s) - - az_ = np.where(top_s > 0, az_ + np.pi, az_) - az_ = np.where(az_ < 0, az_ + 2 * np.pi, az_) - - rg_ = np.sqrt(rx * rx + ry * ry + rz * rz) - el_ = np.arcsin(top_z / rg_) - - return np.rad2deg(az_), np.rad2deg(el_) + return _look_from_ecef( + (pos_x, pos_y, pos_z), + (opos_x, opos_y, opos_z), + lon_rad, + lat_rad, + theta + ) def get_orbit_number(self, utc_time, tbus_style=False, as_float=False): """Calculate orbit number at specified time. @@ -286,19 +362,19 @@ def get_orbit_number(self, utc_time, tbus_style=False, as_float=False): dt = astronomy._days(utc_time - self.orbit_elements.an_time) orbit_period = astronomy._days(self.orbit_elements.an_period) except AttributeError: - pos_epoch, vel_epoch = self.get_position(self.tle.epoch, - normalize=False) + pos_epoch, vel_epoch = self.get_position(self.tle.epoch, normalize=False) if np.abs(pos_epoch[2]) > 1 or not vel_epoch[2] > 0: # Epoch not at ascending node - self.orbit_elements.an_time = self.get_last_an_time( - self.tle.epoch) + self.orbit_elements.an_time = self.get_last_an_time(self.tle.epoch) else: # Epoch at ascending node (z < 1 km) and positive v_z self.orbit_elements.an_time = self.tle.epoch self.orbit_elements.an_period = self.orbit_elements.an_time - \ - self.get_last_an_time(self.orbit_elements.an_time - - np.timedelta64(10, "m")) + self.get_last_an_time(self.orbit_elements.an_time - np.timedelta64(10, "m")) + + logger.debug(f"Orbit reference AN time: {self.orbit_elements.an_time}") + logger.debug(f"Orbit period (days): {astronomy._days(self.orbit_elements.an_period)}") dt = astronomy._days(utc_time - self.orbit_elements.an_time) orbit_period = astronomy._days(self.orbit_elements.an_period) @@ -314,7 +390,16 @@ def get_orbit_number(self, utc_time, tbus_style=False, as_float=False): return orbit - def get_next_passes(self, utc_time, length, lon, lat, alt, tol=0.001, horizon=0): + def get_next_passes( + self, + utc_time: dt.datetime, + length: float, + lon: float, + lat: float, + alt: float, + tol: float = 0.001, + horizon: float = 0 + ) -> list[tuple[dt.datetime, dt.datetime, dt.datetime]]: """Calculate passes for the next hours for a given start time and a given observer. Original by Martin. @@ -331,36 +416,49 @@ def get_next_passes(self, utc_time, length, lon, lat, alt, tol=0.001, horizon=0) """ # every minute - times = utc_time + np.array([dt.timedelta(minutes=minutes) - for minutes in range(length * 60)]) + times = utc_time + np.array([ + dt.timedelta(minutes=minutes) + for minutes in range(round(length * 60)) + ]) elev = self.get_observer_look(times, lon, lat, alt)[1] - horizon zcs = np.where(np.diff(np.sign(elev)))[0] - res = [] - risetime = None - risemins = None + + res: list[tuple[dt.datetime, dt.datetime, dt.datetime]] = [] + risetime: dt.datetime | None = None + risemins: float | None = None + elev_func = partial(self._elevation, utc_time, lon, lat, alt, horizon) elev_inv_func = partial(self._elevation_inv, utc_time, lon, lat, alt, horizon) + for guess in zcs: horizon_mins = _get_root(elev_func, guess, guess + 1.0, tol=tol / 60.0) horizon_time = utc_time + dt.timedelta(minutes=horizon_mins) + if elev[guess] < 0: risetime = horizon_time risemins = horizon_mins else: falltime = horizon_time - fallmins = horizon_mins - if risetime is None: + fallmins: float | None = horizon_mins + + if risetime is None or risemins is None or fallmins is None: continue + int_start = max(0, int(np.floor(risemins))) int_end = min(len(elev), int(np.ceil(fallmins) + 1)) middle = int_start + np.argmax(elev[int_start:int_end]) - highest = utc_time + \ - dt.timedelta(minutes=_get_max_parab( - elev_inv_func, - max(risemins, middle - 1), min(fallmins, middle + 1), - tol=tol / 60.0 - )) - res += [(risetime, falltime, highest)] + + highest = utc_time + dt.timedelta(minutes=_get_max_parab( + elev_inv_func, + max(risemins, middle - 1), + min(fallmins, middle + 1), + tol=tol / 60.0 + )) + + res.append((risetime, falltime, highest)) + risetime = None + risemins = None + return res def _get_time_at_horizon(self, utc_time, obslon, obslat, **kwargs): @@ -488,35 +586,53 @@ def _nprime(time_f): return tcross - - def _elevation(self, utc_time, lon, lat, alt, horizon, minutes): + def _elevation( + self, + utc_time: dt.datetime, + lon: float, + lat: float, + alt: float, + horizon: float, + minutes: float + ) -> float: """Compute the elevation.""" - return self.get_observer_look(utc_time + - dt.timedelta(minutes=np.float64(minutes)), - lon, lat, alt)[1] - horizon - - - def _elevation_inv(self, utc_time, lon, lat, alt, horizon, minutes): + time_f = utc_time + dt.timedelta(minutes=minutes) + _, el = self.get_observer_look(time_f, lon, lat, alt) + return float(el) - horizon + + def _elevation_inv( + self, + utc_time: dt.datetime, + lon: float, + lat: float, + alt: float, + horizon: float, + minutes: float + ) -> float: """Compute the inverse of elevation.""" return -self._elevation(utc_time, lon, lat, alt, horizon, minutes) def _get_root(fun, start, end, tol=0.01): - """Root finding scheme.""" - x_0 = end - x_1 = start - fx_0 = fun(end) - fx_1 = fun(start) + """Find a root of `fun` in [start, end] using Brent's method with tolerance `tol`.""" + x_0 = float(end) + x_1 = float(start) + fx_0 = fun(x_0) + fx_1 = fun(x_1) + + if fx_0 * fx_1 > 0: + raise ValueError("Function values at interval endpoints must have opposite signs.") + + # Swap for better convergence if abs(fx_0) < abs(fx_1): fx_0, fx_1 = fx_1, fx_0 x_0, x_1 = x_1, x_0 - x_n = optimize.brentq(fun, x_0, x_1) - return x_n + return optimize.brentq(fun, x_0, x_1, xtol=tol, rtol=tol) -def _get_max_parab(fun, start, end, tol=0.01): - """Successive parabolic interpolation.""" +def _get_max_parab(fun, start, end, tol=0.01, max_iter=50): + """Find peak of `fun` in [start, end] using parabolic interpolation with fallback.""" a = float(start) c = float(end) b = (a + c) / 2.0 @@ -525,25 +641,35 @@ def _get_max_parab(fun, start, end, tol=0.01): f_b = fun(b) f_c = fun(c) - x = b - with np.errstate(invalid="raise"): - while True: - try: - x = x - 0.5 * (((b - a) ** 2 * (f_b - f_c) - - (b - c) ** 2 * (f_b - f_a)) / - ((b - a) * (f_b - f_c) - (b - c) * (f_b - f_a))) - except FloatingPointError: - return b - if abs(b - x) <= tol: - return x - f_x = fun(x) - # sometimes the estimation diverges... return best guess - if f_x > f_b: - logger.info("Parabolic interpolation did not converge, returning best guess so far.") - return b - - a, b, c = (a + x) / 2.0, x, (x + c) / 2.0 - f_a, f_b, f_c = fun(a), f_x, fun(c) + # Handle flat or symmetric cases + if f_a == f_b == f_c: + return b + + for _ in range(max_iter): + denom = ((b - a) * (f_b - f_c) - (b - c) * (f_b - f_a)) + if denom == 0: + return b # fallback + + try: + x = b - 0.5 * (((b - a)**2 * (f_b - f_c) - (b - c)**2 * (f_b - f_a)) / denom) + except ZeroDivisionError: + return b + + if abs(b - x) <= tol: + return x + + f_x = fun(x) + + # Divergence or invalid result + if f_x > f_b or not (a <= x <= c): + return b + + # Update bracket + a, b, c = (a + x) / 2.0, x, (x + c) / 2.0 + f_a, f_b, f_c = fun(a), f_x, fun(c) + + # Max iterations reached + return b class OrbitElements: @@ -1338,17 +1464,18 @@ def kep2xyz(kep): if __name__ == "__main__": - obs_lon, obs_lat = np.deg2rad((12.4143, 55.9065)) - obs_alt = 0.02 + # Observer's location in degrees + obs_lon, obs_lat = 12.4143, 55.9065 + obs_alt = 0.02 # Altitude in km o = Orbital(satellite="METOP-B") - t_start = dt.datetime.now() - t_stop = t_start + dt.timedelta(minutes=20) + time_step = dt.timedelta(seconds=15) + num_steps = 80 t = t_start - while t < t_stop: - t += dt.timedelta(seconds=15) + print("Time | Azimuth | Elevation | Orbit No.") + for i in range(num_steps): lon, lat, alt = o.get_lonlatalt(t) - lon, lat = np.rad2deg((lon, lat)) az, el = o.get_observer_look(t, obs_lon, obs_lat, obs_alt) ob = o.get_orbit_number(t, tbus_style=True) - print(az, el, ob) + print(f"Time: {t.strftime('%H:%M:%S')} | Az: {az:.2f}° | El: {el:.2f}° | Orbit: {ob}") + t += time_step diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 7ada7d8..940f12b 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -2,12 +2,15 @@ import datetime as dt import unittest +import warnings from unittest import mock +import dask.array as da import numpy as np import pytest +import xarray as xr -from pyorbital import orbital +from pyorbital.orbital import Orbital, get_observer_look eps_deg = 10e-3 @@ -17,22 +20,22 @@ class Test(unittest.TestCase): def test_get_orbit_number(self): """Testing getting the orbitnumber from the TLEs.""" - sat = orbital.Orbital("NPP", - line1="1 37849U 11061A 12017.90990040 " - "-.00000112 00000-0 -32693-4 0 772", - line2="2 37849 98.7026 317.8811 0001845 " - "92.4533 267.6830 14.19582686 11574") + sat = Orbital( + "NPP", + line1="1 37849U 11061A 12017.90990040 -.00000112 00000-0 -32693-4 0 772", + line2="2 37849 98.7026 317.8811 0001845 92.4533 267.6830 14.19582686 11574", + ) dobj = dt.datetime(2012, 1, 18, 8, 4, 19) orbnum = sat.get_orbit_number(dobj) assert orbnum == 1163 def test_sublonlat(self): """Test getting the sub-satellite position.""" - sat = orbital.Orbital("ISS (ZARYA)", - line1="1 25544U 98067A 03097.78853147 " - ".00021906 00000-0 28403-3 0 8652", - line2="2 25544 51.6361 13.7980 0004256 " - "35.6671 59.2566 15.58778559250029") + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) d = dt.datetime(2003, 3, 23, 0, 3, 22) lon, lat, alt = sat.get_lonlatalt(d) expected_lon = -68.199894472013213 @@ -44,11 +47,11 @@ def test_sublonlat(self): def test_observer_look(self): """Test getting the observer look angles.""" - sat = orbital.Orbital("ISS (ZARYA)", - line1="1 25544U 98067A 03097.78853147 " - ".00021906 00000-0 28403-3 0 8652", - line2="2 25544 51.6361 13.7980 0004256 " - "35.6671 59.2566 15.58778559250029") + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) d = dt.datetime(2003, 3, 23, 0, 3, 22) az, el = sat.get_observer_look(d, -84.39733, 33.775867, 0) expected_az = 122.45169655331965 @@ -58,33 +61,33 @@ def test_observer_look(self): def test_orbit_num_an(self): """Test getting orbit number - ascending node.""" - sat = orbital.Orbital("METOP-A", - line1="1 29499U 06044A 11254.96536486 " - ".00000092 00000-0 62081-4 0 5221", - line2="2 29499 98.6804 312.6735 0001758 " - "111.9178 248.2152 14.21501774254058") + sat = Orbital( + "METOP-A", + line1="1 29499U 06044A 11254.96536486 .00000092 00000-0 62081-4 0 5221", + line2="2 29499 98.6804 312.6735 0001758 111.9178 248.2152 14.21501774254058", + ) d = dt.datetime(2011, 9, 14, 5, 30) assert sat.get_orbit_number(d) == 25437 def test_orbit_num_non_an(self): """Test getting orbit number - not ascending node.""" - sat = orbital.Orbital("METOP-A", - line1="1 29499U 06044A 13060.48822809 " - ".00000017 00000-0 27793-4 0 9819", - line2="2 29499 98.6639 121.6164 0001449 " - "71.9056 43.3132 14.21510544330271") + sat = Orbital( + "METOP-A", + line1="1 29499U 06044A 13060.48822809 .00000017 00000-0 27793-4 0 9819", + line2="2 29499 98.6639 121.6164 0001449 71.9056 43.3132 14.21510544330271", + ) dt = np.timedelta64(98, "m") assert sat.get_orbit_number(sat.tle.epoch + dt) == 33028 def test_orbit_num_equator(self): """Test getting orbit numbers when being around equator.""" - sat = orbital.Orbital("SUOMI NPP", - line1="1 37849U 11061A 13061.24611272 " - ".00000048 00000-0 43679-4 0 4334", - line2="2 37849 98.7444 1.0588 0001264 " - "63.8791 102.8546 14.19528338 69643") + sat = Orbital( + "SUOMI NPP", + line1="1 37849U 11061A 13061.24611272 .00000048 00000-0 43679-4 0 4334", + line2="2 37849 98.7444 1.0588 0001264 63.8791 102.8546 14.19528338 69643", + ) t1 = dt.datetime(2013, 3, 2, 22, 2, 25) - t2 = dt.datetime(2013, 3, 2, 22, 2, 26) + t2 = dt.datetime(2013, 3, 2, 22, 3, 00) on1 = sat.get_orbit_number(t1) on2 = sat.get_orbit_number(t2) assert on1 == 6973 @@ -97,25 +100,21 @@ def test_orbit_num_equator(self): def test_get_next_passes_apogee(self): """Regression test #22.""" - line1 = "1 24793U 97020B 18065.48735489 " \ - ".00000075 00000-0 19863-4 0 9994" - line2 = "2 24793 86.3994 209.3241 0002020 " \ - "89.8714 270.2713 14.34246429 90794" + line1 = "1 24793U 97020B 18065.48735489 .00000075 00000-0 19863-4 0 9994" + line2 = "2 24793 86.3994 209.3241 0002020 89.8714 270.2713 14.34246429 90794" - orb = orbital.Orbital("IRIDIUM 7 [+]", line1=line1, line2=line2) + orb = Orbital("IRIDIUM 7 [+]", line1=line1, line2=line2) d = dt.datetime(2018, 3, 7, 3, 30, 15) res = orb.get_next_passes(d, 1, 170.556, -43.368, 0.5, horizon=40) assert abs(res[0][2] - dt.datetime(2018, 3, 7, 3, 48, 13, 178439)) < dt.timedelta(seconds=0.01) def test_get_next_passes_tricky(self): """Check issue #34 for reference.""" - line1 = "1 43125U 18004Q 18251.42128650 " \ - "+.00001666 +00000-0 +73564-4 0 9991" + line1 = "1 43125U 18004Q 18251.42128650 " "+.00001666 +00000-0 +73564-4 0 9991" - line2 = "2 43125 097.5269 314.3317 0010735 "\ - "157.6344 202.5362 15.23132245036381" + line2 = "2 43125 097.5269 314.3317 0010735 " "157.6344 202.5362 15.23132245036381" - orb = orbital.Orbital("LEMUR-2-BROWNCOW", line1=line1, line2=line2) + orb = Orbital("LEMUR-2-BROWNCOW", line1=line1, line2=line2) d = dt.datetime(2018, 9, 8) res = orb.get_next_passes(d, 72, -8.174163, 51.953319, 0.05, horizon=5) @@ -130,9 +129,9 @@ def test_get_next_passes_issue_22(self): line1 = "1 28654U 05018A 21083.16603416 .00000102 00000-0 79268-4 0 9999" line2 = "2 28654 99.0035 147.6583 0014816 159.4931 200.6838 14.12591533816498" - orb = orbital.Orbital("NOAA 18", line1=line1, line2=line2) + orb = Orbital("NOAA 18", line1=line1, line2=line2) t = dt.datetime(2021, 3, 9, 22) - next_passes = orb.get_next_passes(t, 1, -15.6335, 27.762, 0.) + next_passes = orb.get_next_passes(t, 1, -15.6335, 27.762, 0.0) rise, fall, max_elevation = next_passes[0] assert rise < max_elevation < fall @@ -140,17 +139,18 @@ def test_get_next_passes_issue_22(self): def test_utc2local(self, get_lonlatalt): """Test converting UTC to local time.""" get_lonlatalt.return_value = -45, None, None - sat = orbital.Orbital("METOP-A", - line1="1 29499U 06044A 13060.48822809 " - ".00000017 00000-0 27793-4 0 9819", - line2="2 29499 98.6639 121.6164 0001449 " - "71.9056 43.3132 14.21510544330271") + sat = Orbital( + "METOP-A", + line1="1 29499U 06044A 13060.48822809 .00000017 00000-0 27793-4 0 9819", + line2="2 29499 98.6639 121.6164 0001449 71.9056 43.3132 14.21510544330271", + ) assert sat.utc2local(dt.datetime(2009, 7, 1, 12)) == dt.datetime(2009, 7, 1, 9) @mock.patch("pyorbital.orbital.Orbital.utc2local") @mock.patch("pyorbital.orbital.Orbital.get_orbit_number") def test_get_equatorial_crossing_time(self, get_orbit_number, utc2local): """Test get the equatorial crossing time.""" + def get_orbit_number_patched(utc_time, **kwargs): utc_time = np.datetime64(utc_time) diff = (utc_time - np.datetime64("2009-07-01 12:38:12")) / np.timedelta64(7200, "s") @@ -158,31 +158,386 @@ def get_orbit_number_patched(utc_time, **kwargs): get_orbit_number.side_effect = get_orbit_number_patched utc2local.return_value = "local_time" - sat = orbital.Orbital("METOP-A", - line1="1 29499U 06044A 13060.48822809 " - ".00000017 00000-0 27793-4 0 9819", - line2="2 29499 98.6639 121.6164 0001449 " - "71.9056 43.3132 14.21510544330271") + sat = Orbital( + "METOP-A", + line1="1 29499U 06044A 13060.48822809 .00000017 00000-0 27793-4 0 9819", + line2="2 29499 98.6639 121.6164 0001449 71.9056 43.3132 14.21510544330271", + ) # Ascending node - res = sat.get_equatorial_crossing_time(tstart=dt.datetime(2009, 7, 1, 12), - tend=dt.datetime(2009, 7, 1, 13)) + res = sat.get_equatorial_crossing_time(tstart=dt.datetime(2009, 7, 1, 12), tend=dt.datetime(2009, 7, 1, 13)) exp = dt.datetime(2009, 7, 1, 12, 38, 12) assert res - exp < dt.timedelta(seconds=0.01) # Descending node - res = sat.get_equatorial_crossing_time(tstart=dt.datetime(2009, 7, 1, 12), - tend=dt.datetime(2009, 7, 1, 14, 0), - node="descending") + res = sat.get_equatorial_crossing_time( + tstart=dt.datetime(2009, 7, 1, 12), + tend=dt.datetime(2009, 7, 1, 14, 0), + node="descending", + ) exp = dt.datetime(2009, 7, 1, 13, 38, 12) assert res - exp < dt.timedelta(seconds=0.01) # Conversion to local time - res = sat.get_equatorial_crossing_time(tstart=dt.datetime(2009, 7, 1, 12), - tend=dt.datetime(2009, 7, 1, 14), - local_time=True) + res = sat.get_equatorial_crossing_time( + tstart=dt.datetime(2009, 7, 1, 12), + tend=dt.datetime(2009, 7, 1, 14), + local_time=True, + ) assert res == "local_time" + def test_get_next_passes_basic(self): + """Test basic pass detection with known TLE and observer.""" + line1 = "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652" + line2 = "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029" + sat = Orbital("ISS (ZARYA)", line1=line1, line2=line2) + + observer_lon = -84.39733 + observer_lat = 33.775867 + observer_alt = 0.0 + start_time = dt.datetime(2003, 3, 23, 0, 0) + + passes = sat.get_next_passes(start_time, length=2, lon=observer_lon, lat=observer_lat, alt=observer_alt) + assert len(passes) > 0 + for rise, fall, peak in passes: + assert rise < peak < fall + + def test_get_next_passes_no_passes(self): + """Test when no passes occur within the time window.""" + line1 = "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652" + line2 = "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029" + sat = Orbital("ISS (ZARYA)", line1=line1, line2=line2) + + observer_lon = 0.0 + observer_lat = -90.0 # South Pole, unlikely for ISS + observer_alt = 0.0 + start_time = dt.datetime(2003, 3, 23, 0, 0) + + passes = sat.get_next_passes(start_time, length=1, lon=observer_lon, lat=observer_lat, alt=observer_alt) + assert passes == [] + + def test_get_next_passes_with_horizon(self): + """Test pass detection with a non-zero horizon elevation.""" + line1 = "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652" + line2 = "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029" + sat = Orbital("ISS (ZARYA)", line1=line1, line2=line2) + + observer_lon = -84.39733 + observer_lat = 33.775867 + observer_alt = 0.0 + start_time = dt.datetime(2003, 3, 23, 0, 0) + + passes = sat.get_next_passes( + start_time, + length=2, + lon=observer_lon, + lat=observer_lat, + alt=observer_alt, + horizon=20, + ) + for rise, fall, peak in passes: + assert rise < peak < fall + + def test_get_next_passes_tolerance_effect(self): + """Test that changing tolerance affects timing precision.""" + line1 = "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652" + line2 = "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029" + sat = Orbital("ISS (ZARYA)", line1=line1, line2=line2) + + observer_lon = -84.39733 + observer_lat = 33.775867 + observer_alt = 0.0 + start_time = dt.datetime(2003, 3, 23, 0, 0) + + passes_lo_tol = sat.get_next_passes( + start_time, + length=2, + lon=observer_lon, + lat=observer_lat, + alt=observer_alt, + tol=0.001, + ) + passes_hi_tol = sat.get_next_passes( + start_time, + length=2, + lon=observer_lon, + lat=observer_lat, + alt=observer_alt, + tol=1.0, + ) + + assert len(passes_lo_tol) == len(passes_hi_tol) + for (r1, f1, p1), (r2, f2, p2) in zip(passes_lo_tol, passes_hi_tol): + assert abs((r1 - r2).total_seconds()) < 60 # Should be close, but not identical + + def test_observer_look_radian_bug(self): + """Detect bug caused by passing observer coordinates in radians instead of degrees.""" + line1 = "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652" + line2 = "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029" + sat = Orbital("ISS (ZARYA)", line1=line1, line2=line2) + + time = dt.datetime(2003, 3, 23, 0, 3, 22) + obs_lon_deg, obs_lat_deg = -84.39733, 33.775867 + obs_lon_rad, obs_lat_rad = np.deg2rad(obs_lon_deg), np.deg2rad(obs_lat_deg) + + # Correct usage: degrees + az_deg, el_deg = sat.get_observer_look(time, obs_lon_deg, obs_lat_deg, 0) + + # Incorrect usage: radians + az_rad, el_rad = sat.get_observer_look(time, obs_lon_rad, obs_lat_rad, 0) + + assert abs(az_deg - az_rad) > 10, "Azimuth unexpectedly similar—possible bug hidden" + assert abs(el_deg - el_rad) > 10, "Elevation unexpectedly similar—possible bug hidden" + + def test_get_next_passes_radian_bug(self): + """Detect bug caused by passing observer coordinates in radians instead of degrees.""" + line1 = "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652" + line2 = "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029" + sat = Orbital("ISS (ZARYA)", line1=line1, line2=line2) + + start_time = dt.datetime(2003, 3, 23, 0, 0) + observer_lon_deg, observer_lat_deg = -84.39733, 33.775867 + observer_lon_rad, observer_lat_rad = np.deg2rad(observer_lon_deg), np.deg2rad(observer_lat_deg) + + # Correct usage: degrees + passes_deg = sat.get_next_passes(start_time, length=2, lon=observer_lon_deg, lat=observer_lat_deg, alt=0) + + # Incorrect usage: radians + passes_rad = sat.get_next_passes(start_time, length=2, lon=observer_lon_rad, lat=observer_lat_rad, alt=0) + + assert len(passes_deg) > 0, "Expected passes not found with correct input" + assert len(passes_rad) == 0 or any( + abs((p1[0] - p2[0]).total_seconds()) > 600 for p1, p2 in zip(passes_deg, passes_rad) + ), "Radian input unexpectedly produced similar results—possible bug hidden" + + def test_find_aos_basic(self): + """Test basic AOS detection for a known satellite and location.""" + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) + utc_time = dt.datetime(2003, 3, 23, 0, 0, 0) + aos = sat.find_aos(utc_time, lon=-84.39733, lat=33.775867, alt=0) + assert aos is not None, "AOS should be detected" + assert aos > utc_time, "AOS must be after the input time" + + def test_find_aol_basic(self): + """Test basic AOL detection for a known satellite and location.""" + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) + utc_time = dt.datetime(2003, 3, 23, 0, 0, 0) + aol = sat.find_aol(utc_time, lon=-84.39733, lat=33.775867, alt=0) + assert aol is not None, "AOL should be detected" + assert aol > utc_time, "AOL must be after the input time" + + def test_find_aos_no_pass(self): + """Test AOS detection when no pass occurs within search window.""" + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) + # Use a location far from satellite's inclination + utc_time = dt.datetime(2003, 3, 23, 0, 0, 0) + aos = sat.find_aos(utc_time, lon=0, lat=89.0, alt=0, horizon=10) + assert aos is None, "No AOS should be found near the poles for ISS" + + def test_find_aos_edge_horizon(self): + """Test AOS detection with a high horizon angle.""" + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) + utc_time = dt.datetime(2003, 3, 23, 0, 0, 0) + aos = sat.find_aos(utc_time, lon=-84.39733, lat=33.775867, alt=0, horizon=30) + assert aos is None or aos > utc_time, "AOS may not occur with high horizon angle" + + +@pytest.mark.parametrize( + ("utc_time", "lon", "lat", "alt", "horizon", "is_aos", "expect_result"), + [ + # Basic pass over Atlanta + (dt.datetime(2003, 3, 23, 0, 0, 0), -84.39733, 33.775867, 0, 0, True, True), + # High elevation requirement (still detectable) + (dt.datetime(2003, 3, 23, 0, 0, 0), -84.39733, 33.775867, 0, 30, True, True), + # Near polar region (unlikely pass for ISS) + (dt.datetime(2003, 3, 23, 0, 0, 0), 0, 89.0, 0, 0, True, False), + # AOL detection instead of AOS + (dt.datetime(2003, 3, 23, 0, 0, 0), -84.39733, 33.775867, 0, 0, False, True), + # Satellite below horizon at start time (should still detect future AOS) + (dt.datetime(2003, 3, 23, 0, 0, 0), -84.39733, 33.775867, 0, 10, True, True), + # Observer at sea level vs. high altitude (compare AOS detectability) + (dt.datetime(2003, 3, 23, 0, 0, 0), -84.39733, 33.775867, 3000, 10, True, False), + # Horizon angle just below satellite max elevation (should fail) + (dt.datetime(2003, 3, 23, 0, 0, 0), -84.39733, 33.775867, 0, 89, True, False), + # Negative altitude (below sea level, e.g., Dead Sea) + (dt.datetime(2003, 3, 23, 0, 0, 0), 35.5, 31.5, -430, 10, True, True), + # Equator, zero horizon — should always detect AOS if satellite is visible + (dt.datetime(2003, 3, 23, 0, 0, 0), 0, 0, 0, 0, True, True), + ], +) +def test_crossing_detection(utc_time, lon, lat, alt, horizon, is_aos, expect_result): + """Test horizon crossing detection under various observer conditions.""" + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) + + result = ( + sat.find_aos(utc_time, lon, lat, alt, horizon) if is_aos else sat.find_aol(utc_time, lon, lat, alt, horizon) + ) + + if expect_result: + assert result is not None, "Expected crossing, got None" + assert result > utc_time, "Crossing must be after input time" + else: + assert result is None, "Expected no crossing, but got one" + + +@pytest.mark.parametrize( + ("utc_time", "lon", "lat", "alt", "horizon", "is_aos", "expect_result"), + [ + # Extreme case: ISS cannot reach 89° latitude with 80° elevation mask + (dt.datetime(2003, 3, 23, 0, 0, 0), 0, 89.0, 0, 80, True, False), + ], +) +def test_extreme_no_aos(utc_time, lon, lat, alt, horizon, is_aos, expect_result): + """Test that no AOS is detected under extreme elevation and location constraints.""" + sat = Orbital( + "ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ) + + result = ( + sat.find_aos(utc_time, lon, lat, alt, horizon) if is_aos else sat.find_aol(utc_time, lon, lat, alt, horizon) + ) + + if expect_result: + assert result is not None, "Expected crossing, got None" + assert result > utc_time, "Crossing must be after input time" + else: + assert result is None, "Expected no crossing, but got one" + + +@pytest.mark.parametrize( + ("tle", "test_time", "expected_sign"), + [ + # ISS: should be near equator, Z < 0 before ascending node + ( + ( + "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ), + dt.datetime(2003, 3, 23, 0, 0, 0), + 1, # Expect Z velocity > 0 at ascending node + ), + # METOP-A: polar orbit, ascending node should be near equator + ( + ( + "1 29499U 06044A 13060.48822809 .00000017 00000-0 27793-4 0 9819", + "2 29499 98.6639 121.6164 0001449 71.9056 43.3132 14.21510544330271", + ), + dt.datetime(2013, 3, 1, 12, 0, 0), + 1, + ), + # SUOMI NPP: test near epoch + ( + ( + "1 37849U 11061A 13061.24611272 .00000048 00000-0 43679-4 0 4334", + "2 37849 98.7444 1.0588 0001264 63.8791 102.8546 14.19528338 69643", + ), + dt.datetime(2013, 3, 2, 6, 0, 0), + 1, + ), + ], +) +def test_get_last_an_time_velocity_sign(tle, test_time, expected_sign): + """Test that get_last_an_time returns a time with positive Z velocity (ascending node).""" + sat = Orbital("TestSat", line1=tle[0], line2=tle[1]) + an_time = sat.get_last_an_time(test_time) + _, vel = sat.get_position(an_time, normalize=False) + assert vel[2] * expected_sign > 0, f"Expected Z velocity sign {expected_sign}, got {vel[2]}" + + +@pytest.mark.parametrize( + ("tle", "test_time"), + [ + # ISS + ( + ( + "1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + "2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029", + ), + dt.datetime(2003, 3, 23, 0, 0, 0), + ), + # METOP-A + ( + ( + "1 29499U 06044A 13060.48822809 .00000017 00000-0 27793-4 0 9819", + "2 29499 98.6639 121.6164 0001449 71.9056 43.3132 14.21510544330271", + ), + dt.datetime(2013, 3, 1, 12, 0, 0), + ), + ], +) +def test_get_last_an_time_z_velocity_positive(tle, test_time): + """Verify that Z velocity is positive at the computed ascending node time.""" + sat = Orbital("TestSat", line1=tle[0], line2=tle[1]) + an_time = sat.get_last_an_time(test_time) + _, vel = sat.get_position(an_time, normalize=False) + assert vel[2] > 0, f"Expected positive Z velocity, got {vel[2]}" + + +@pytest.mark.parametrize( + ("tle", "test_time"), + [ + # SUOMI NPP + ( + ( + "1 37849U 11061A 13061.24611272 .00000048 00000-0 43679-4 0 4334", + "2 37849 98.7444 1.0588 0001264 63.8791 102.8546 14.19528338 69643", + ), + dt.datetime(2013, 3, 2, 6, 0, 0), + ), + ], +) +def test_get_last_an_time_period_consistency(tle, test_time): + """Ensure the time between consecutive ascending nodes matches the orbital period.""" + sat = Orbital("TestSat", line1=tle[0], line2=tle[1]) + an1 = sat.get_last_an_time(test_time) + an2 = sat.get_last_an_time(an1 - np.timedelta64(10, "m")) + period = (an1 - an2).astype("timedelta64[s]").item().total_seconds() / 60.0 + expected_period = 1440.0 / sat.tle.mean_motion # minutes + assert abs(period - expected_period) < 1.0, f"Period mismatch: got {period}, expected {expected_period}" + + +@pytest.mark.parametrize( + ("tle", "test_time"), + [ + # METOP-A near node + ( + ( + "1 29499U 06044A 13060.48822809 .00000017 00000-0 27793-4 0 9819", + "2 29499 98.6639 121.6164 0001449 71.9056 43.3132 14.21510544330271", + ), + dt.datetime(2013, 3, 1, 5, 30, 0), + ), + ], +) +def test_get_last_an_time_near_node(tle, test_time): + """Confirm ascending node time is reasonably close to the input time.""" + sat = Orbital("TestSat", line1=tle[0], line2=tle[1]) + an_time = sat.get_last_an_time(test_time) + delta = abs((np.datetime64(test_time) - an_time).astype("timedelta64[s]").item().total_seconds()) + assert delta < 7200, f"AN time too far from test time: {delta} seconds" + class TestGetObserverLook(unittest.TestCase): """Test the get_observer_look function.""" @@ -190,91 +545,85 @@ class TestGetObserverLook(unittest.TestCase): def setUp(self): """Set up the test environment.""" self.t = dt.datetime(2018, 1, 1, 0, 0, 0) - self.sat_lon = np.array([[-89.5, -89.4, -89.5, -89.4], - [-89.3, -89.2, -89.3, -89.2]]) - self.sat_lat = np.array([[45.5, 45.4, 45.5, 45.4], - [45.3, 40.2, 45.3, 40.2]]) + self.sat_lon = np.array([[-89.5, -89.4, -89.5, -89.4], [-89.3, -89.2, -89.3, -89.2]]) + self.sat_lat = np.array([[45.5, 45.4, 45.5, 45.4], [45.3, 40.2, 45.3, 40.2]]) self.sat_alt = 35786 * np.ones((2, 4)) - self.lon = np.array([[-85.5, -85.4, -89.5, -99.4], - [-85.3, -89.2, -89.3, -79.2]]) - self.lat = np.array([[40.5, 40.4, 65.5, 45.4], - [40.3, 40.2, 25.3, 40.2]]) + self.lon = np.array([[-85.5, -85.4, -89.5, -99.4], [-85.3, -89.2, -89.3, -79.2]]) + self.lat = np.array([[40.5, 40.4, 65.5, 45.4], [40.3, 40.2, 25.3, 40.2]]) self.alt = np.zeros((2, 4)) - self.exp_azi = np.array([[331.00275902, 330.95954165, 180, 86.435411], - [330.91642994, 180, 0, 273.232073]]) - self.exp_elev = np.array([[83.18070976, 83.17788976, 66.548467, 81.735221], - [83.17507167, 90, 66.559906, 81.010018]]) + self.exp_azi = np.array( + [ + [331.00275902, 330.95954165, 180, 86.435411], + [330.91642994, 180, 0, 273.232073], + ] + ) + self.exp_elev = np.array( + [ + [83.18070976, 83.17788976, 66.548467, 81.735221], + [83.17507167, 90, 66.559906, 81.010018], + ] + ) def test_basic_numpy(self): """Test with numpy array inputs.""" - from pyorbital import orbital - azi, elev = orbital.get_observer_look(self.sat_lon, self.sat_lat, - self.sat_alt, self.t, - self.lon, self.lat, self.alt) + azi, elev = get_observer_look( + self.sat_lon, + self.sat_lat, + self.sat_alt, + self.t, + self.lon, + self.lat, + self.alt, + ) np.testing.assert_allclose(azi, self.exp_azi) np.testing.assert_allclose(elev, self.exp_elev) def test_basic_dask(self): """Test with dask array inputs.""" - import dask.array as da - - from pyorbital import orbital sat_lon = da.from_array(self.sat_lon, chunks=2) sat_lat = da.from_array(self.sat_lat, chunks=2) sat_alt = da.from_array(self.sat_alt, chunks=2) lon = da.from_array(self.lon, chunks=2) lat = da.from_array(self.lat, chunks=2) alt = da.from_array(self.alt, chunks=2) - azi, elev = orbital.get_observer_look(sat_lon, sat_lat, - sat_alt, self.t, - lon, lat, alt) + azi, elev = get_observer_look(sat_lon, sat_lat, sat_alt, self.t, lon, lat, alt) np.testing.assert_allclose(azi.compute(), self.exp_azi) np.testing.assert_allclose(elev.compute(), self.exp_elev) def test_xarray_with_numpy(self): """Test with xarray DataArray with numpy array as inputs.""" - import xarray as xr - - from pyorbital import orbital def _xarr_conv(input_array): return xr.DataArray(input_array) + sat_lon = _xarr_conv(self.sat_lon) sat_lat = _xarr_conv(self.sat_lat) sat_alt = _xarr_conv(self.sat_alt) lon = _xarr_conv(self.lon) lat = _xarr_conv(self.lat) alt = _xarr_conv(self.alt) - azi, elev = orbital.get_observer_look(sat_lon, sat_lat, - sat_alt, self.t, - lon, lat, alt) + azi, elev = get_observer_look(sat_lon, sat_lat, sat_alt, self.t, lon, lat, alt) np.testing.assert_allclose(azi.data, self.exp_azi) np.testing.assert_allclose(elev.data, self.exp_elev) def test_xarray_with_dask(self): """Test with xarray DataArray with dask array as inputs.""" - import dask.array as da - import xarray as xr - - from pyorbital import orbital def _xarr_conv(input_array): return xr.DataArray(da.from_array(input_array, chunks=2)) + sat_lon = _xarr_conv(self.sat_lon) sat_lat = _xarr_conv(self.sat_lat) sat_alt = _xarr_conv(self.sat_alt) lon = _xarr_conv(self.lon) lat = _xarr_conv(self.lat) alt = _xarr_conv(self.alt) - azi, elev = orbital.get_observer_look(sat_lon, sat_lat, - sat_alt, self.t, - lon, lat, alt) + azi, elev = get_observer_look(sat_lon, sat_lat, sat_alt, self.t, lon, lat, alt) np.testing.assert_allclose(azi.data.compute(), self.exp_azi) - np.testing.assert_allclose(elev.data.compute(), self.exp_elev) + np.testing.assert_allclose(elev, self.exp_elev) def test_scalar(self): """Test with scalar inputs.""" - from pyorbital.orbital import get_observer_look (azi, elev) = get_observer_look(0, 0, 30_000_000, self.t, 0, 0, 0) np.testing.assert_allclose(elev, 90) @@ -304,74 +653,65 @@ def setUp(self): def test_basic_numpy(self): """Test with numpy array inputs.""" - from pyorbital import orbital - azi, elev = orbital.get_observer_look(self.sat_lon, self.sat_lat, - self.sat_alt, self.t, - self.lon, self.lat, self.alt) + azi, elev = get_observer_look( + self.sat_lon, + self.sat_lat, + self.sat_alt, + self.t, + self.lon, + self.lat, + self.alt, + ) assert np.sum(np.isnan(azi)) == 0 assert not np.isnan(azi).any() np.testing.assert_allclose(elev, self.exp_elev) def test_basic_dask(self): """Test with dask array inputs.""" - import dask.array as da - - from pyorbital import orbital sat_lon = da.from_array(self.sat_lon, chunks=2) sat_lat = da.from_array(self.sat_lat, chunks=2) sat_alt = da.from_array(self.sat_alt, chunks=2) lon = da.from_array(self.lon, chunks=2) lat = da.from_array(self.lat, chunks=2) alt = da.from_array(self.alt, chunks=2) - azi, elev = orbital.get_observer_look(sat_lon, sat_lat, - sat_alt, self.t, - lon, lat, alt) + azi, elev = get_observer_look(sat_lon, sat_lat, sat_alt, self.t, lon, lat, alt) assert np.sum(np.isnan(azi)) == 0 assert not np.isnan(azi).any() np.testing.assert_allclose(elev.compute(), self.exp_elev) def test_xarray_with_numpy(self): """Test with xarray DataArray with numpy array as inputs.""" - import xarray as xr - - from pyorbital import orbital def _xarr_conv(input_array): return xr.DataArray(input_array) + sat_lon = _xarr_conv(self.sat_lon) sat_lat = _xarr_conv(self.sat_lat) sat_alt = _xarr_conv(self.sat_alt) lon = _xarr_conv(self.lon) lat = _xarr_conv(self.lat) alt = _xarr_conv(self.alt) - azi, elev = orbital.get_observer_look(sat_lon, sat_lat, - sat_alt, self.t, - lon, lat, alt) + azi, elev = get_observer_look(sat_lon, sat_lat, sat_alt, self.t, lon, lat, alt) assert np.sum(np.isnan(azi)) == 0 assert not np.isnan(azi).any() np.testing.assert_allclose(elev.data, self.exp_elev) def test_xarray_with_dask(self): """Test with xarray DataArray with dask array as inputs.""" - import dask.array as da - import xarray as xr - - from pyorbital import orbital def _xarr_conv(input_array): return xr.DataArray(da.from_array(input_array, chunks=2)) + sat_lon = _xarr_conv(self.sat_lon) sat_lat = _xarr_conv(self.sat_lat) sat_alt = _xarr_conv(self.sat_alt) lon = _xarr_conv(self.lon) lat = _xarr_conv(self.lat) alt = _xarr_conv(self.alt) - azi, elev = orbital.get_observer_look(sat_lon, sat_lat, - sat_alt, self.t, - lon, lat, alt) + azi, elev = get_observer_look(sat_lon, sat_lat, sat_alt, self.t, lon, lat, alt) assert np.sum(np.isnan(azi)) == 0 assert not np.isnan(azi).any() - np.testing.assert_allclose(elev.data.compute(), self.exp_elev) + np.testing.assert_allclose(elev, self.exp_elev) class TestRegressions(unittest.TestCase): @@ -379,45 +719,50 @@ class TestRegressions(unittest.TestCase): def test_63(self): """Check that no runtimewarning is raised, #63.""" - import warnings - - from pyorbital.orbital import Orbital warnings.filterwarnings("error") - orb = Orbital("Suomi-NPP", - line1="1 37849U 11061A 19292.84582509 .00000011 00000-0 25668-4 0 9997", - line2="2 37849 98.7092 229.3263 0000715 98.5313 290.6262 14.19554485413345") + orb = Orbital( + "Suomi-NPP", + line1="1 37849U 11061A 19292.84582509 .00000011 00000-0 25668-4 0 9997", + line2="2 37849 98.7092 229.3263 0000715 98.5313 290.6262 14.19554485413345", + ) orb.get_next_passes(dt.datetime(2019, 10, 21, 16, 0, 0), 12, 123.29736, -13.93763, 0) warnings.filterwarnings("default") -@pytest.mark.parametrize("dtime", - [dt.datetime(2024, 6, 25, 11, 0, 18), - dt.datetime(2024, 6, 25, 11, 5, 0, 0, dt.timezone.utc), - np.datetime64("2024-06-25T11:10:00.000000") - ] - ) +@pytest.mark.parametrize( + "dtime", + [ + dt.datetime(2024, 6, 25, 11, 0, 18), + dt.datetime(2024, 6, 25, 11, 5, 0, 0, dt.timezone.utc), + np.datetime64("2024-06-25T11:10:00.000000"), + ], +) def test_get_last_an_time_scalar_input(dtime): """Test getting the time of the last ascending node - input time is a scalar.""" - from pyorbital.orbital import Orbital - orb = Orbital("NOAA-20", - line1="1 43013U 17073A 24176.73674251 .00000000 00000+0 11066-3 0 00014", - line2="2 43013 98.7060 114.5340 0001454 139.3958 190.7541 14.19599847341971") + orb = Orbital( + "NOAA-20", + line1="1 43013U 17073A 24176.73674251 .00000000 00000+0 11066-3 0 00014", + line2="2 43013 98.7060 114.5340 0001454 139.3958 190.7541 14.19599847341971", + ) expected = np.datetime64("2024-06-25T10:44:18.234375") result = orb.get_last_an_time(dtime) assert abs(expected - result) < np.timedelta64(1, "s") -@pytest.mark.parametrize("dtime", - [dt.datetime(2024, 6, 25, 11, 5, 0, 0, dt.timezone(dt.timedelta(hours=1))), - ] - ) +@pytest.mark.parametrize( + "dtime", + [ + dt.datetime(2024, 6, 25, 11, 5, 0, 0, dt.timezone(dt.timedelta(hours=1))), + ], +) def test_get_last_an_time_wrong_input(dtime): """Test getting the time of the last ascending node - wrong input.""" - from pyorbital.orbital import Orbital - orb = Orbital("NOAA-20", - line1="1 43013U 17073A 24176.73674251 .00000000 00000+0 11066-3 0 00014", - line2="2 43013 98.7060 114.5340 0001454 139.3958 190.7541 14.19599847341971") + orb = Orbital( + "NOAA-20", + line1="1 43013U 17073A 24176.73674251 .00000000 00000+0 11066-3 0 00014", + line2="2 43013 98.7060 114.5340 0001454 139.3958 190.7541 14.19599847341971", + ) expected = "UTC time expected! Parsing a timezone aware datetime object requires it to be UTC!" with pytest.raises(ValueError, match=expected): diff --git a/pyorbital/tests/test_utils.py b/pyorbital/tests/test_utils.py new file mode 100644 index 0000000..ecc6772 --- /dev/null +++ b/pyorbital/tests/test_utils.py @@ -0,0 +1,218 @@ +"""Test suite for utils in pyorbital.orbital.""" + +import numpy as np +import pytest + +from pyorbital.orbital import _get_max_parab, _get_root, compute_azimuth_elevation, ecef_to_topocentric + + +@pytest.mark.parametrize( + ("fun", "a", "b", "expected"), + [ + (lambda x: x**2 - 4, 1, 3, 2.0), # Basic root test + (lambda x: x**2 - 4, -3, -1, -2.0), # Negative root test + (lambda x: np.sin(x), 3, 4, np.pi), # High precision near pi + ], +) +def test_get_root_parametrized(fun, a, b, expected): + """Test root finding with various functions and intervals.""" + root = _get_root(fun, a, b, tol=1e-6) + assert np.isclose(root, expected, atol=1e-3) + + +def test_get_root_unbracketed(): + """Test that root finding fails when no root is bracketed.""" + + def fun(x): + return x**2 + 1 + + with pytest.raises(ValueError, match=".*opposite signs.*"): + _get_root(fun, -1, 1) + + +@pytest.mark.parametrize( + ("fun", "a", "b", "expected", "atol"), + [ + (lambda x: -(x**2) + 4 * x, 1, 3, 2.0, 1e-3), # Basic peak test + (lambda x: 5, 0, 10, 5.0, 1.0), # Flat function + (lambda x: -np.cos(x), 0, np.pi, np.pi / 2, 1e-2), # Cosine peak + ], +) +def test_get_max_parab_parametrized(fun, a, b, expected, atol): + """Test peak finding with various functions and intervals.""" + peak = _get_max_parab(fun, a, b) + assert np.isclose(peak, expected, atol=atol) + + +def test_get_max_parab_divergence(): + """Test fallback behavior for a rising exponential with no peak.""" + + def fun(x): + return np.exp(x) + + peak = _get_max_parab(fun, 0, 1) + assert 0 <= peak <= 1 + + +def test_ecef_to_topocentric_scalar(): + """Test topocentric conversion with scalar inputs.""" + rx, ry, rz = 1.0, 2.0, 3.0 + lat, lon, theta = np.deg2rad(45), np.deg2rad(60), np.deg2rad(30) + top_s, top_e, top_z = ecef_to_topocentric(rx, ry, rz, lat, lon, theta) + + assert isinstance(top_s, float) + assert isinstance(top_e, float) + assert isinstance(top_z, float) + + +def test_ecef_to_topocentric_array(): + """Test topocentric conversion with array inputs.""" + rx = np.array([1.0, 2.0]) + ry = np.array([2.0, 3.0]) + rz = np.array([3.0, 4.0]) + lat = np.deg2rad(np.array([45.0, 46.0])) + lon = np.deg2rad(np.array([60.0, 61.0])) + theta = np.deg2rad(np.array([30.0, 31.0])) + + top_s, top_e, top_z = ecef_to_topocentric(rx, ry, rz, lat, lon, theta) + assert top_s.shape == (2,) + assert top_e.shape == (2,) + assert top_z.shape == (2,) + + +def test_compute_azimuth_elevation_scalar(): + """Test azimuth and elevation computation with scalar inputs.""" + top_s, top_e, top_z = 1.0, 1.0, 1.0 + rg = np.sqrt(3) + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + + assert isinstance(az, float) + assert isinstance(el, float) + assert 0 <= az <= 360 + assert -90 <= el <= 90 + + +@pytest.mark.parametrize( + ("top_s", "top_e", "top_z"), + [ + ( + np.array([1.0, 0.0]), + np.array([0.0, 1.0]), + np.array([1.0, 1.0]), + ), # Basic array test + ], +) +def test_compute_azimuth_elevation_array(top_s, top_e, top_z): + """Test azimuth and elevation computation with array inputs.""" + rg = np.sqrt(top_s**2 + top_e**2 + top_z**2) + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + assert az.shape == (2,) + assert el.shape == (2,) + assert np.all((az >= 0) & (az <= 360)) + assert np.all((el >= -90) & (el <= 90)) + + +def test_elevation_clipping(): + """Test elevation clipping when top_z slightly exceeds rg.""" + top_s = np.array([0.0]) + top_e = np.array([0.0]) + top_z = np.array([1.0000001]) + rg = np.array([1.0]) + + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + assert np.isclose(el, 90.0) + + +def test_zero_distance(): + """Test behavior when observer and satellite are at the same location.""" + rx = ry = rz = 0.0 + rg = 1e-10 # Avoid division by zero + az, el = compute_azimuth_elevation(rx, ry, rz, rg) + assert np.isclose(el, 0.0) or np.isclose(el, 90.0) + + +def test_negative_altitude(): + """Validate azimuth/elevation output when observer is below sea level.""" + rx, ry, rz = 1.0, 1.0, 1.0 + lat, lon, theta = np.deg2rad(0), np.deg2rad(0), np.deg2rad(0) + top_s, top_e, top_z = ecef_to_topocentric(rx, ry, rz, lat, lon, theta) + rg = np.sqrt(rx**2 + ry**2 + rz**2) + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + assert 0 <= az <= 360 + assert -90 <= el <= 90 + + +def test_polar_observer(): + """Check azimuth/elevation computation near the geographic poles.""" + lat = np.deg2rad(89.999) + lon = np.deg2rad(0) + theta = np.deg2rad(0) + rx, ry, rz = 1.0, 0.0, 0.0 + top_s, top_e, top_z = ecef_to_topocentric(rx, ry, rz, lat, lon, theta) + rg = np.sqrt(rx**2 + ry**2 + rz**2) + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + assert 0 <= az <= 360 + assert -90 <= el <= 90 + + +def test_large_array_input(): + """Stress test azimuth/elevation computation with large input arrays.""" + n = 10000 + rx = np.ones(n) + ry = np.ones(n) + rz = np.ones(n) + lat = np.deg2rad(np.full(n, 45.0)) + lon = np.deg2rad(np.full(n, 60.0)) + theta = np.deg2rad(np.full(n, 30.0)) + top_s, top_e, top_z = ecef_to_topocentric(rx, ry, rz, lat, lon, theta) + rg = np.sqrt(rx**2 + ry**2 + rz**2) + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + assert az.shape == (n,) + assert el.shape == (n,) + + +def test_nan_inf_input(): + """Ensure azimuth/elevation handles NaN and Inf inputs gracefully.""" + rx = np.array([np.nan, np.inf]) + ry = np.array([1.0, 1.0]) + rz = np.array([1.0, 1.0]) + rg = np.sqrt(rx**2 + ry**2 + rz**2) + az, el = compute_azimuth_elevation(rx, ry, rz, rg) + assert np.isnan(az[0]) + assert np.isfinite(az[1]) + assert np.isnan(el[0]) + assert np.isfinite(el[1]) + + +def test_empty_array_input(): + """Ensure azimuth/elevation handles zero-length arrays gracefully.""" + top_s = np.array([]) + top_e = np.array([]) + top_z = np.array([]) + rg = np.array([]) + + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + assert az.size == 0 + assert el.size == 0 + + +@pytest.mark.parametrize( + ("lat_deg", "lon_deg"), + [ + (90.0, 0.0), # North Pole + (-90.0, 0.0), # South Pole + (0.0, 180.0), # International Date Line East + (0.0, -180.0), # International Date Line West + ], +) +def test_extreme_lat_lon(lat_deg, lon_deg): + """Test topocentric conversion at extreme latitude/longitude values.""" + rx, ry, rz = 1.0, 1.0, 1.0 + lat = np.deg2rad(lat_deg) + lon = np.deg2rad(lon_deg) + theta = 0.0 + top_s, top_e, top_z = ecef_to_topocentric(rx, ry, rz, lat, lon, theta) + rg = np.sqrt(rx**2 + ry**2 + rz**2) + az, el = compute_azimuth_elevation(top_s, top_e, top_z, rg) + assert 0 <= az <= 360 + assert -90 <= el <= 90