Bugfix use of datetime.datetime objects in call to get_last_an_time

pyorbital/orbital.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2011-2023 Pyorbital developers
+# Copyright (c) 2011-2024 Pyorbital developers
 
 # Author(s):
 
@@ -27,6 +27,7 @@
 import logging
 import warnings
 from datetime import datetime, timedelta
+import pytz
 
 import numpy as np
 from scipy import optimize
@@ -86,11 +87,14 @@
 
 
 class OrbitalError(Exception):
+    """Custom exception for the Orbital class."""
+
     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)
@@ -142,7 +146,6 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt):
 
 
 class Orbital(object):
-
     """Class for orbital computations.
 
     The *satellite* parameter is the name of the satellite to work on and is
@@ -151,6 +154,7 @@ class Orbital(object):
     """
 
     def __init__(self, satellite, tle_file=None, line1=None, line2=None):
+        """Initialize the class."""
         satellite = satellite.upper()
         self.satellite_name = satellite
         self.tle = tlefile.read(satellite, tle_file=tle_file,
@@ -159,16 +163,14 @@ def __init__(self, satellite, tle_file=None, line1=None, line2=None):
         self._sgdp4 = _SGDP4(self.orbit_elements)
 
     def __str__(self):
+        """Print the Orbital object state."""
         return self.satellite_name + " " + str(self.tle)
 
     def get_last_an_time(self, utc_time):
-        """Calculate time of last ascending node relative to the
-        specified time
-        """
-
+        """Calculate time of last ascending node relative to the specified time."""
         # Propagate backwards to ascending node
         dt = np.timedelta64(10, 'm')
-        t_old = utc_time
+        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)
@@ -236,23 +238,26 @@ def get_lonlatalt(self, utc_time):
         return np.rad2deg(lon), np.rad2deg(lat), alt
 
     def find_aos(self, utc_time, lon, lat):
+        """Find AOS."""
         pass
 
     def find_aol(self, utc_time, lon, lat):
+        """Find AOL."""
         pass
 
     def get_observer_look(self, utc_time, lon, lat, alt):
         """Calculate observers look angle to a satellite.
-        http://celestrak.com/columns/v02n02/
+
+        See 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)
-        """
 
+        """
         utc_time = dt2np(utc_time)
         (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = self.get_position(
             utc_time, normalize=False)
@@ -330,8 +335,7 @@ 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):
-        """Calculate passes for the next hours for a given start time and a
-        given observer.
+        """Calculate passes for the next hours for a given start time and a given observer.
 
         Original by Martin.
 
@@ -344,8 +348,8 @@ def get_next_passes(self, utc_time, length, lon, lat, alt, tol=0.001, horizon=0)
         :horizon: the elevation of horizon to compute risetime and falltime.
 
         :return: [(rise-time, fall-time, max-elevation-time), ...]
-        """
 
+        """
         def elevation(minutes):
             """Compute the elevation."""
             return self.get_observer_look(utc_time +
@@ -358,7 +362,7 @@ def elevation_inv(minutes):
             return -elevation(minutes)
 
         def get_root(fun, start, end, tol=0.01):
-            """Root finding scheme"""
+            """Root finding scheme."""
             x_0 = end
             x_1 = start
             fx_0 = fun(end)
@@ -432,15 +436,18 @@ def get_max_parab(fun, start, end, tol=0.01):
         return res
 
     def _get_time_at_horizon(self, utc_time, obslon, obslat, **kwargs):
-        """Get the time closest in time to *utc_time* when the
-        satellite is at the horizon relative to the position of an observer on
-        ground (altitude = 0)
+        """Determine when the satellite is at the horizon relative to an observer on ground.
+
+        Get the time closest in time to *utc_time* when the satellite is at the
+        horizon relative to the position of an observer on ground (altitude =
+        0).
 
         Note: This is considered deprecated and it's functionality is currently
         replaced by 'get_next_passes'.
+
         """
         warnings.warn("_get_time_at_horizon is replaced with get_next_passes",
-                      DeprecationWarning)
+                      DeprecationWarning, stacklevel=2)
         if "precision" in kwargs:
             precision = kwargs['precision']
         else:
@@ -518,7 +525,7 @@ def get_equatorial_crossing_time(self, tstart, tend, node='ascending', local_tim
             return None
         elif n_end - n_start > 1:
             warnings.warn('Multiple revolutions between start and end time. Computing crossing '
-                          'time for the last revolution in that interval.')
+                          'time for the last revolution in that interval.', stacklevel=2)
 
         # Let n'(t) = n(t) - offset. Determine offset so that n'(tstart) < 0 and n'(tend) > 0 and
         # n'(tcross) = 0.
@@ -555,11 +562,10 @@ def _nprime(time_f):
 
 
 class OrbitElements(object):
-
-    """Class holding the orbital elements.
-    """
+    """Class holding the orbital elements."""
 
     def __init__(self, tle):
+        """Initialize the class."""
         self.epoch = tle.epoch
         self.excentricity = tle.excentricity
         self.inclination = np.deg2rad(tle.inclination)
@@ -609,9 +615,7 @@ def __init__(self, tle):
 
 
 class _SGDP4(object):
-
-    """Class for the SGDP4 computations.
-    """
+    """Class for the SGDP4 computations."""
 
     def __init__(self, orbit_elements):
         self.mode = None
@@ -917,7 +921,25 @@ def propagate(self, utc_time):
         return kep
 
 
+def _get_tz_unaware_utctime(utc_time):
+    """Return timzone unaware datetime object.
+
+    The input *utc_time* is either a timezone unaware object assumed to be in
+    UTC, or a timezone aware datetime object in UTC.
+    """
+    if isinstance(utc_time, datetime):
+        if utc_time.tzinfo and utc_time.tzinfo != pytz.utc:
+            raise ValueError("UTC time expected! Parsing a timezone aware datetime object requires it to be UTC!")
+        return utc_time.replace(tzinfo=None)
+
+    return utc_time
+
+
 def kep2xyz(kep):
+    """Keppler to cartesian coordinates conversion.
+
+    (Not sure what 'kep' actually refers to, just guessing! FIXME!)
+    """
     sinT = np.sin(kep['theta'])
     cosT = np.cos(kep['theta'])
     sinI = np.sin(kep['eqinc'])