Issues predicting satellite positions

[a-tangible-human-skull]

I am attempting to simulate the imagery generated by a wide FOV telescope imager, so I can do some overlays over my observations. As such, I am generating the positions of stars and satellites at an epoch, determining their relative positions to my telescope in RA and DEC, and then projecting them into camera coordinate space to get their apparent pixel position given my telescope's orientation.

Stars work great! But for some reason, I get very bizarre behavior when I try to add satellites to the test patterns.

I'm sure I'm messing up something about the methodology, or misinterpreting how epochs are handled when pulling data from the propagated satellite catalog, but for the life of me I've yet been unable to get satellites to align well with my observations. While no doubt this stems in part with the challenges of getting the epoch of the imagery correct, the behavior of satellites in my synthetic imagery makes me think I'm interpreting some input/output incorrectly.

Here's a reduced snippet of my code:

import os
from math import cos, sin, radians

from PIL import Image, ImageDraw
from skyfield.api import load, Star, wgs84, N, W, Topos
from skyfield.data import hipparcos
from pytz import timezone

from datetime import datetime, timedelta

# -------------------------------------------------
# Define camera paremeters, imaged part of sky, etc

camera_state = {'lat': 32.2319,
                'lon': 110.9501,

                'width_pix': 1280,
                'height_pix': 720,

                'fov_deg': 44.7,

                'ra_deg': 279.864,
                'dec_deg': 32.809,
                'roll_deg': -90.0
                }

magnitude_thresh = 6.0

camera_fov_rad = radians(camera_state['fov_deg'])

# Doing attitude in RA/DEC/ROLL for alignment with star tracking for orientation ID
camera_att = {'ra': radians(camera_state['ra_deg']),
              'dec': radians(camera_state['dec_deg']),
              'roll': radians(camera_state['roll_deg'])}

camera_center = (camera_state['width_pix'] / 2,
                 camera_state['height_pix'] / 2)

ifov = (camera_state['width_pix'] / camera_fov_rad,
        camera_state['width_pix'] / camera_fov_rad)  # Square pixels

# Epoch for observations
local_date = datetime(year=2023,
                      month=8,
                      day=24,
                      hour=20,
                      minute=1,
                      second=8,
                      microsecond=0)

loc_timezone = timezone('US/Pacific')

# Large full sky 3les
catalog_path = r"G:\Projects\Video\Misc\astro\20230824\3le.txt"

# ------------------------------------------------
# Parameters for outputs
output_path = r'G:\Projects\Video\Misc\astro\20230824\test_output'
timestep_count = 5

# -------------------------------------------------
# Load ephemeris

planets = load('de421.bsp')
earth = planets['earth']
observor_loc_w_earth = earth + wgs84.latlon(camera_state['lat'] * N, camera_state['lon'] * W)
observor_loc_no_earth = wgs84.latlon(camera_state['lat'] * N, camera_state['lon'] * W)  # Have this here in case I want to try the alt method
observor_alt_method = Topos(latitude=str(camera_state['lat']) + ' ' + 'N', longitude=str(camera_state['lon']) + ' ' + 'W')

# Stars
with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)
filtered_stars_df = stars[stars['magnitude'] <= magnitude_thresh]
print('Filtered stars by magnitude %s from %s to %s' % (magnitude_thresh, len(stars), len(filtered_stars_df)))

# Satellites
sats = None
if os.path.exists(catalog_path):
    sats = load.tle_file(catalog_path)


# ---------------------------------------------------
# Evaluate if the observed elements are in the FOV
# Terminology is star based, but should hold for satellites as well...

def predict_fov_position(star_ra_deg, star_dec_deg, camera_att_in, ifov_in, camera_state_in, camera_center_in):
    # Conversions
    s_ra = radians(star_ra_deg)
    s_dec = radians(star_dec_deg)

    # Coordinates
    stdX = cos(s_dec) * sin(s_ra - camera_att_in['ra']) / \
           (cos(camera_att_in['dec']) * cos(s_dec) * cos(s_ra - camera_att_in['ra']) + sin(camera_att_in['dec']) * sin(s_dec))

    stdY = (sin(camera_att_in['dec']) * cos(s_dec) * cos(s_ra - camera_att_in['ra']) - cos(camera_att_in['dec']) * sin(s_dec)) / \
           (cos(camera_att_in['dec']) * cos(s_dec) * cos(s_ra - camera_att_in['ra']) + sin(camera_att_in['dec']) * sin(s_dec))

    # Calculate the pixel coordinates of the star in the image
    starPixelX = cos(camera_att_in['roll']) * ifov_in[0] * stdX + sin(camera_att_in['roll']) * ifov_in[1] * stdY + camera_center_in[0]
    starPixelY = -sin(camera_att_in['roll']) * ifov_in[0] * stdX + cos(camera_att_in['roll']) * ifov_in[1] * stdY + camera_center_in[1]
    starPixelX = camera_state_in['width_pix'] - starPixelX

    return starPixelX, starPixelY


# ----------------------------------------------------
stars_in_fov = []  # Never going to change
framecount = 0
# Loop through some number of time steps so satellite motion can be observed
for timestep in range(0, timestep_count, 1):

    sats_in_fov = []
    framecount = framecount + 1

    print('Timestep %s' % (timestep,))

    # Offset epoch
    time_delta = timedelta(seconds=float(timestep))
    curr_epoch = loc_timezone.localize(local_date + time_delta)
    ts = load.timescale()
    timescale = ts.from_datetime(curr_epoch)

    # ----------------------------------------------------
    if framecount == 1:  # Stars are the same for each frame given our lack of attitude shift here
        # Loop through star catalog
        for i in filtered_stars_df.index:
            current_star = Star.from_dataframe(filtered_stars_df.loc[i])

            astrometric = observor_loc_w_earth.at(timescale).observe(current_star)
            ra, dec, distance = astrometric.radec()

            ra_deg = ra._degrees
            dec_deg = dec.degrees

            x_out, y_out = predict_fov_position(ra_deg, dec_deg, camera_att, ifov, camera_state, camera_center)

            if (camera_state['width_pix'] > x_out) and (camera_state['height_pix'] > y_out) and (x_out >= 0) and (y_out >= 0):
                stars_in_fov.append(
                    {
                        'cat_id': i,
                        'location': [x_out, y_out],
                        'magnitude': filtered_stars_df.loc[i].magnitude
                    }
                )

        print("Number of stars in FOV of %s is: %s" % (camera_state['fov_deg'], len(stars_in_fov)))

    # ----------------------------------------------------
    # Loop through sat catalog
    if sats is not None:
        for sat in sats:

            # Method 1
            orbit = (sat - observor_alt_method).at(timescale)
            ra, dec, distance = orbit.radec()
            alt, az, distance = orbit.altaz()

            # Method 2
            '''
            difference = sat - observor_loc_no_earth
            topocentric = difference.at(timescale)
            ra, dec, distance = topocentric.radec() # Using RA/Dec for alignment with star method
            alt, az, distance = topocentric.altaz()
            '''
            ra_deg = ra._degrees
            dec_deg = dec.degrees

            # Remove satellites not within a certain distance
            height = wgs84.height_of(sat.at(timescale)).km
            if (height > 1000):
                continue

            # Remove if below the Horizon?
            if alt.degrees < 0.0:
                continue

            # Remove satellites that aren't sunlit
            sunlit = sat.at(timescale).is_sunlit(planets)  # This also doesn't seem to be working right...
            if not sunlit:
                continue

            x_out, y_out = predict_fov_position(ra_deg, dec_deg, camera_att, ifov, camera_state, camera_center)

            if (camera_state['width_pix'] > x_out) and (camera_state['height_pix'] > y_out) and (x_out >= 0) and (y_out >= 0):
                sats_in_fov.append(
                    {
                        'name': sat.name,
                        'cat_id': sat.model.satnum,
                        'location': [x_out, y_out]
                    }
                )

    print("Number of sats in FOV of %s is: %s" % (camera_state['fov_deg'], len(sats_in_fov)))

    # -----------------------
    # Generate imagery

    extendedPoint_size = 2.0

    imagery = Image.new("RGB", (camera_state['width_pix'], camera_state['height_pix']), "black")
    draw = ImageDraw.Draw(imagery)

    # Loop through stars
    for star in stars_in_fov:
        fill = 'white'
        if (star['magnitude'] < 3.0):  # Make brighter stars bigger for eyeballing
            useEllipse = True
        else:
            useEllipse = False

        if not useEllipse:
            draw.point((star['location'][0], star['location'][1]), fill=fill)
        else:
            draw.ellipse((star['location'][0] - extendedPoint_size, star['location'][1] - extendedPoint_size,
                          star['location'][0] + extendedPoint_size, star['location'][1] + extendedPoint_size),
                         fill=fill, outline=fill)

    # Loop through sats
    for sat in sats_in_fov:
        draw.ellipse((sat['location'][0] - extendedPoint_size, sat['location'][1] - extendedPoint_size,
                      sat['location'][0] + extendedPoint_size, sat['location'][1] + extendedPoint_size),
                     fill="#ffff33", outline="#ffff33")

    # Save imagery
    out_path = r"%s\ra%s_dec%s_roll%s_fov%s_mag%s_epoch_%s.png" % (output_path,
                                                                   camera_state['ra_deg'],
                                                                   camera_state['dec_deg'],
                                                                   camera_state['roll_deg'],
                                                                   camera_state['fov_deg'],
                                                                   magnitude_thresh,
                                                                   timescale.astimezone(timezone('US/Pacific')).strftime("%m/%d/%Y, %H:%M:%S.%f").replace(
                                                                       ':', '-').replace('/', '-'))

    imagery.save(out_path)

Absurd amounts of satellites show up as flagged as both being sunlit and in my camera FOV. Thousands. No amount of filtering seems to get rid of them either; altitude, sunlit, etc. Additionally, I get satellites whose altitude should be below the horizon supposedly having RA/Decs that would put them in my FOV, which... doesn't make any sense.

I'm also almost positive my satellite altitude check is just straight up wrong, since it's reporting no satellites below 400 km in my FOV which... seems unlikely given my wide FOV.

There's got to be something really stupidly obvious that I'm missing, but I'm not seeing it. I appreciate any insight.

[brandon-rhodes]

Do you get exactly the same results with your # Method 2? I can't tell from your description whether the two methods behave the same.

[a-tangible-human-skull]

They produce identical results; the output images are binary identical.

I had implemented both approaches in an attempt to debug this behavior, and left them in just case I was missing a step in either methodology.

edit: as a general demonstration that the stellar approach here seems to work, here's an image from my WFOV imagery and a compared image. The differences between the two should stem from issues in approximating the camera FOV, which I can solve by using the ESA tetra3 library on a bunch of subsequent stills to produce a better FOV estimate. There's also certain stars missing from either due to band differences between my camera and Hipparcos. Also, I'm only capturing at a 1/30 sec integration time so a lot of stuff falls behind noise.

[brandon-rhodes]

Three quick Saturday thoughts, in case any are helpful:

1. Ah, I'm happy to hear the two approaches produce the same image; otherwise we would have other mysteries to solve!
2. Alas, I can't see any correlation between the stars in the image (where I see the W of Cassiopeia, I think?) and your synthetic star image. I would be tempted to keep working on the projection routine until I saw a clear W of exactly the same shape; otherwise, the problem could be the projection.
3. How many satellites are in your satellite file? Your field of view is pretty big, and all other things being equal, if you're showing, say, 1/30th of the sky, you'd expect to see 1/30th of the sky's satellites. Except — except for the fact that you're up on one side of the Earth, which places very few satellites overhead (except maybe out in the ring of geostationary satellites, which I think you're filtering out?) because at any given moment most low-earth-orbiting satellites are "behind" you behind the horizon. But a question arises: what if you are accidentally showing all the satellites behind you, down beneath the Earth? That would be a huge number of satellites. Maybe you could start spot-checking a few satellites and see; my thought is that maybe they are behind you somehow.

[a-tangible-human-skull]

Ah blast, I think you're right; passing the generated image into Astrometry.net as a sanity check returns an attitude that is inverted compared to the input. I'm pretty sure I had double checked this as working a few weeks ago, so I must have broken something in the interim... As you say, that would cause the phenomenology I'm seeing here, especially since I am using a full space-track catalog. I'll keep poking. Thanks!

[brandon-rhodes]

I should admit up front that I don't have a good story for why if alt.degrees < 0.0: wouldn't filter out most satellites, regardless of any problems with the projection, since that if condition doesn't depend on the behavior of the projection. But without the files involved I can't try out the code myself, so I'm left trying to guess. I'll let you know if I think of anything else!

[a-tangible-human-skull]

OK, so good news: everything seems to be working now. Stars are oriented correctly relative to my observations, and I am even seeing satellites in my emulated scene that appear in my observations. Which is fantastic! I'm not entirely sure what I meaningfully changed that made such a difference beyond adjusting the emulated camera settings to match closer to my observations.

I am seeing some slight deltas in satellite positions relative to the stars that seem large relative to the anticipated catalog errors (which, AFAIK, should be only a few kilometers off, which would be basically sub-pixel with this camera's giant IFOV). Some of that is tied to epoch and the fact that I don't think I have my RA/DEC exactly right, but I feel like there's probably some other, additional element I'm missing.

Here's a comparison:

Red dots are the skyfield stars. Yellow dots are skyfield satellites. I overlaid a transparency of my observations, and colored the satellite observations I can make out teal. Note in particular the fairly consistent vertical bias on satellites, perpendicular to the satellite motion vector.

Updated code:
unit_test_sats.txt

Here's a 3le file limited to the objects in the FOV:
3le_191.txt

The bright object in the center of the FOV is the RESURS DK-1, and STARLINK-5849 is the faintly visible in the lower right in my imagery.

I'll have to go through some other imagery I took to see if this bias is consistent beyond this limited data set, but is there anything obvious in my handling of the skyfield satellites in particular that would cause these kinds of deltas? The level of delta here would correspond to a massive catalog error, which I doubt.

[brandon-rhodes]

I'll have to go through some other imagery I took to see if this bias is consistent beyond this limited data set, but is there anything obvious in my handling of the skyfield satellites in particular that would cause these kinds of deltas?

No, I don't see anything wrong with your satellite handling; it seems to match what's in Skyfield's documentation.

I'm a bit confused about why most of the red Skyfield stars don't match up with stars in your photograph. Why isn't Skyfield showing the W of Cassiopeia? Do its stars fall below the threshold you're setting for magnitude?