Refactor orbital tracking with vectorized methods, bug fixes and improved utilities

[JaskRendix]

PR improves the orbital tracking module and related utilities.

• fixed a bug in the main script where observer longitude and latitude were being converted to radians before being passed to get_observer_look, which already performs this conversion internally, coordinates are now passed in degrees as expected, and conversion is handled inside the method

• refactored the get_observer_look function: split the logic into two helper functions: ecef_to_topocentric and compute_azimuth_elevation, this isolates coordinate transformations and makes the azimuth/elevation computation easier to test and reuse, pytest coverage added

• refactored _get_root and _get_max_parab: _get_root now includes a check for valid root bracketing and uses xtol and rtol for better convergence control. _get_max_parab now includes a fallback for flat functions, a maximum iteration limit, and improved numerical stability to avoid divergence and handle edge cases.

• cleaned up test_orbital: ran isort and black for formatting and added new pytest coverage for all newly introduced vectorized methods and utilities

• implemented find_aos and find_aol methods for detecting Acquisition of Signal and Loss of Signal events based on observer location and elevation mask; includes robust horizon-crossing logic and full pytest coverage

• improved get_last_an_time logic by refactoring shared node detection into a private _find_last_node_time method, which uses direct geometric analysis of the satellite’s Z-position and velocity to identify node crossings; this approach is distinct from orbit-number-based equatorial crossing methods and supports both ascending and descending detection, added get_last_dn_time to retrieve the last descending node, addressing feature request Would be nice to have a get_last_dn_time to get the last descending node #95; includes parameterized pytest coverage for velocity sign, node classification, and temporal accuracy

• Closes Would be nice to have a get_last_dn_time to get the last descending node #95

• Tests added

• Fully documented

[djhoese]

Nice. At a high level I have two big questions:

1. Is it possible to make the original methods work vectorized without needed separate "vectorized" specific methods?
2. Would using something like numpy.vectorize help? https://numpy.org/doc/stable/reference/generated/numpy.vectorize.html. I think in this case the outputs would always be numpy arrays but I think that's OK.

One other concern is the use of np.array to convert to array types and isinstance checks on numpy arrays. It'd be nice if these functions were compatible with dask arrays which likely (hopefully) just treating things like numpy arrays without specifically checking. There are also functions like https://numpy.org/doc/2.2/reference/generated/numpy.asanyarray.html#numpy.asanyarray which may help with being a little more flexible.

In the long run I suppose it'd be nice to have the low-level functions support multiple values instead of having to iterate over everything every time. I think I tried doing that once but decided it was too hard.

Thoughts?

[codecov]

Codecov Report

❌ Patch coverage is 95.70136% with 19 lines in your changes missing coverage. Please review.
✅ Project coverage is 91.06%. Comparing base (c2cbf0c) to head (58cbc63).
⚠️ Report is 2 commits behind head on main.

Files with missing lines	Patch %	Lines
pyorbital/orbital.py	87.68%	17 Missing ⚠️
pyorbital/tests/test_orbital.py	98.87%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #206      +/-   ##
==========================================
+ Coverage   90.46%   91.06%   +0.60%     
==========================================
  Files          19       20       +1     
  Lines        3040     3313     +273     
==========================================
+ Hits         2750     3017     +267     
- Misses        290      296       +6     

Flag	Coverage Δ
unittests	91.06% <95.70%> (+0.60%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[coveralls]

coverage: 91.052% (+0.6%) from 90.431%
when pulling e68719d on JaskRendix:vectorized
into cd410a9 on pytroll:main.

[mraspaud]

I'm pretty sure get_lonlatalt at least already support arrays as input, what is get_lonlatalt_vectorized offering more?

[JaskRendix]

@mraspaud absolutely nothing, I'm going to remove these

edit:

I’ve removed all of the _vectorized wrapper methods. While they offered an array-friendly API, they didn't provide true performance. I'm convinced that the SGP4 calculation is the real performance constraint and fixing it there will automatically make all public methods correctly array-aware.

[djhoese]

Thanks for responding so quickly to the reviews! Do these current changes change the existing return type to any of these methods?

Have you (could you) run mypy on your the files you've changed? I don't think we currently run any type checkers on pyorbital automatically (CI or pre-commit) so it'd be good to know that mypy agrees.

[JaskRendix]

@djhoese I usually run mypy --strict, but it throws a ton of errors. Do you use a specific mypy command? If not, I can just strip out the typing so it matches the rest of the code, no big deal.

I checked and:
get_lonlatalt still returns (float,float,float) or a tuple of np.ndarray's if vectorized
get_last_an_time still returns a single time object (np.datetime64)
_elevation still returns a single float
_elevation_inv still returns a single float
_get_root still returns a single float
_get_max_parab still returns a single float
both get_observer_look still returns a tuple of two floats

_find_single_crossing, find_aos, find_aol, ecef_to_topocentric and compute_azimuth_elevation are new

I'm going to run mypy less strict soon

[djhoese]

Typically I just go with mypy <pkgdir>/, but pyorbital is well behind other pytroll packages on linting/typing standards so we don't have anything existing setup.

[JaskRendix]

@djhoese yeah, I created 2 typehints, but I have fixed both

pyorbital/orbital.py:458: error: Argument 1 to "__call__" of "_UFunc_Nin1_Nout1" has incompatible type "Any | None"; expected "complex | str | bytes | generic[Any]" [arg-type]
pyorbital/orbital.py:464: error: Value of type variable "SupportsRichComparisonT" of "max" cannot be "signedinteger[_32Bit | _64Bit] | Any | None" [type-var]

now back to 5

(.venv) giorgio@giorgio-Aspire-V3-572G:~/pyorbital$ mypy pyorbital
pyorbital/orbital.py:42: error: Incompatible types in assignment (expression has type "None", variable has type Module)  [assignment]
pyorbital/orbital.py:49: error: Incompatible types in assignment (expression has type "None", variable has type Module)  [assignment]
pyorbital/tests/test_astronomy.py:38: error: Cannot assign to a type  [misc]
pyorbital/tests/test_astronomy.py:38: error: Incompatible types in assignment (expression has type "None", variable has type "type[DataArray]")  [assignment]
pyorbital/geoloc_example.py:29: error: Skipping analyzing "mpl_toolkits.basemap": module is installed, but missing library stubs or py.typed marker  [import-untyped]
pyorbital/geoloc_example.py:29: note: See https://mypy.readthedocs.io/en/stable/running_mypy.html#missing-imports
Found 5 errors in 3 files (checked 20 source files)

[djhoese]

Ok those errors are from ImportError hackery. I'm OK ignore them for now and in the future we can add an inline type ignore comment.

@mraspaud @pnuu do either of you have time to run Satpy's unit tests against this PR and make sure it passes so we're not surprised later?

[pnuu]

I'll have a go running the Satpy tests.

[mraspaud]

I also have tests I can run against other libs

[djhoese]

Unrelated: @JaskRendix this project is planned to be relicensed from GPL to Apache version 2. Once we organize some of the other dependencies we'll ask all past contributor's if they're OK with this in a GitHub issue. This is just a heads up.

[pnuu]

With Satpy I only get an unrelated failure in downloading TLEs from Celestrak, most likely our network is being throttled or something. So all good from Satpy's perspective!

[pnuu]

Pytroll-collectors and Trollflow2 are fine, but I got a failure in Pytroll-Schedule in https://github.com/pytroll/pytroll-schedule/blob/main/trollsched/tests/test_schedule.py#L221 , the expected rise time and resulting value differ by few milliseconds. Might be close enough 😅 This is most likely coming from the bug you fixed in the previous PR. I'll fix this in Pytroll-Schedule, we might not need sub-second accuracy in scheduling...

rt1 = datetime.datetime(2018, 11, 28, 10, 53, 42, 79483)
rise_times = [datetime.datetime(2018, 11, 28, 10, 53, 42, 66585),
                       datetime.datetime(2018, 11, 28, 12, 34, 44, 662042)]

[djhoese]

Ref: pytroll/pytroll-schedule#110

[djhoese]

Now that pytroll-schedule is cleaned up, do we have any other issues with this PR? Merge?

[mraspaud]

I haven't tested it yet, I'll try to get to it tonight

[djhoese]

Ah sorry, I misread your previous comment to mean you had tested those things.

[mraspaud]

Thanks a lot for this PR, I appreciate the huge refactoring effort!

I have some comments and questions inline.

[mraspaud]

why remove the compute call here and not on the line above?

[JaskRendix]

Because azi is still a dask‑backed DataArray, so you must trigger computation, while elev is already a plain NumPy array coming out of the math path and has nothing to compute.

[mraspaud]

Would it be possible to use the Orbital's get_observer_look to refactor some of this function.

[JaskRendix]

I pulled the shared ECEF > topocentric math into a small helper and had both functions call it, so the standalone version keeps its lon/lat/alt semantics and Orbital.get_observer_look keeps its TLE‑based propagation, but the duplicated geometry is gone.

[mraspaud]

I understand this is more secure to ensure the loop ever ends, thanks a lot!

One question I have is how was the default value of 50 chosen?

[JaskRendix]

50‑iteration cap wasn’t chosen for any deep numerical reason, it’s a safe upper bound that guarantees the loop can’t run forever. Parabolic interpolation settles in well under 10–15 steps, so 50 is above any realistic need. We could drop it to something like 20 without changing the behavior in normal cases.

[mraspaud]

why this except clause if we already test above that denom is 0?

[JaskRendix]

Because the denom == 0 check only catches the exact 0.0 case, while the division can still blow up on things like -0.0 or values so close to zero that floating‑point quirks trigger a ZeroDivisionError. The except is just a safety net for those edge cases.