Issue Calculating Accurate View Zenith Angles on Terra Satellite Overpasses

Hello!

I am working on processing the View Zenith Angles (VZA) of the MODIS instrument onboard the Terra satellite. I'm targeting specific latitude and longitude coordinates at predetermined times. My goal is to determine the closest overpasses with a VZA <= ~60°, to align with the maximum VZA of the MODIS instrument, within a 24/48 hour timeframe. Basically I would like to find the satellite overpass with the lowest VZA within the last 24 hours, meaning the satellite that was as close to directly overhead as possible.

Currently my code gives me output values of elevation (with which I compute VZA = 90° - abs(elevation)) that are very small (always close to 0) which makes me question how those satellite overpasses are determined. This would mean that MODIS did not capture any of these points, which is quite unlikely in practice. I wonder whether there is an issue with how I'm using the get_next_passes and get_observer_look functions, or perhaps I am misunderstanding some of the settings for filtering passes based on elevation. I tried filtering elevation between different values but it does not change the values. When I set the horizon to a specific value, the observations generally cluster around those values.

Does someone know how the Pyorbital overpass method works, or what I might be doing wrong?

Here is my code with a few coordinates set for a single time, which should be replicable:

!pip install pyorbital

from pyorbital.orbital import Orbital
from datetime import datetime, timedelta
import pandas as pd

# Defining latitudes and longitudes of points of interest
coordinates = [
    (34.634998, -105.23556),
    (33.558334, -81.791664),
    (30.283333, -85.783333),
    (61.766666, -150.60001), 
    (65.066666, -146.16667)
]

date = datetime.fromisoformat("2020-12-23T00:00:00")
date = date.replace(hour=14, minute=13)
date = date - timedelta(days=1)

# Initialize the orbital object for a satellite using its TLE data
line1 = '1 25994U 99068A   20357.07230622  .00000069  00000+0  25285-4 0  9990'
line2 = '2 25994  98.1966  68.5759 0001201  98.2528  18.8483 14.57114596117631'
orb = Orbital("Terra", line1=line1, line2=line2)

# Initialize a list to collect results for each coordinate
results = []

# Iterate over each pair of latitude and longitude
for idx, (latitude, longitude) in enumerate(coordinates):
    print(f"Checking passes for coordinates: Latitude {latitude}, Longitude {longitude}")
    pass_times = orb.get_next_passes(date, 24, longitude, latitude, alt=0, tol=0.001, horizon=0)

    if not pass_times:
        print(f"No visible passes for coordinates: Latitude {latitude}, Longitude {longitude}")
    else:
        print(f"Found {len(pass_times)} passes")

    max_elevation = -90  # Initialize to a very low value to ensure any valid pass is considered
    best_pass = None

    for pass_time in pass_times:
        rise_time, max_time, set_time = pass_time
        azimuth, elevation = orb.get_observer_look(max_time, longitude, latitude, 0)
        print(f"At {max_time}, Azimuth {azimuth}, Elevation {elevation}")

        # Ensure elevation is strictly between 30 and 90 degrees
        if 30 <= elevation <= 90:
            if elevation > max_elevation:  # Update best pass if this pass has a higher elevation
                max_elevation = elevation
                best_pass = (azimuth, elevation, max_time)

    if best_pass:
        results.append({
            "location": idx,
            "latitude": latitude,
            "longitude": longitude,
            "azimuth": best_pass[0],
            "elevation": best_pass[1],
            "passtime": best_pass[2]
        })

# Check the content of results list
if results:
    print("Results were found and will be displayed.")
    rel_data = pd.DataFrame(results)
    print(rel_data)
else:
    print("No results were found, no passes met the criteria.")

Any help or guidance would be greatly appreciated!

The passlist returns a list of tuples where the first and second elements are the start and end of the local pass relative to your "station", thus you are looking at the two horizons, the rise and set of the satellite.

You would need to pass the third element of the 3-tuple to get_observer_look

In your example above this is wrong! rise_time, max_time, set_time = pass_time

It should be: rise_time, set_time, max_time = pass_time

I have done the same mistake myself before. I think the order of the elements are not intuitive, perhaps something we should fix?

Example:

In [29]: aqua = Orbital('EOS-Aqua')
In [30]: passlist = aqua.get_next_passes(now, 6, 16, 59, 0, horizon=0)
In [31]: passlist
Out[31]: 
[(datetime.datetime(2024, 5, 7, 23, 53, 34, 204465),
  datetime.datetime(2024, 5, 8, 0, 5, 7, 636094),
  datetime.datetime(2024, 5, 7, 23, 59, 22, 348750)),
 (datetime.datetime(2024, 5, 8, 1, 30, 55, 637588),
  datetime.datetime(2024, 5, 8, 1, 44, 53, 264153),
  datetime.datetime(2024, 5, 8, 1, 37, 56, 314023))]

In [32]: aqua.get_observer_look(passlist[0][2], 16, 59, 0)
Out[32]: (82.34616159635466, 13.514754629087017)

pytroll / pyorbital

Issue Calculating Accurate View Zenith Angles on Terra Satellite Overpasses #151