Topos errors when calculating visible EarthSatellite

[Kahn]

I feel like this is going to be my fault but here we go. I'm working on a python 2.7.x codebase and found AttributeError: 'Topos' object has no attribute 'position' when trying to calculate visible EarthSatellites. I've tested again with a clean Python 3.7 environment with the same result though =\

The use case is mapping Iridium NEXT orbits against known positions and times.

If anyone has any pointers or can see where I have flubbed that would be awesome!

Reproduction

Platform: Windows 10 python-3.7.1 jplephem-2.9 numpy-1.16.1 sgp4-1.4 skyfield-1.10

Code

from skyfield.api import Topos, load

ts = load.timescale()
t = ts.utc(2019, 2, 17, 11, 18, 7)
stations_url = 'http://celestrak.com/NORAD/elements/iridium-NEXT.txt'
satellites = load.tle(stations_url)
print(satellites['IRIDIUM 161'])
sat_pos_geocentric = satellites['IRIDIUM 161'].at(t)
print(sat_pos_geocentric.position.km)
subpoint = sat_pos_geocentric.subpoint()
print('Latitude:', subpoint.latitude)
print('Longitude:', subpoint.longitude)
print('Elevation (m):', int(subpoint.elevation.m))
bluffton = Topos('40.8939 N', '83.8917 W')
difference = sat_pos_geocentric - bluffton
print(difference)

Output

$ python skyfield-test.py
EarthSatellite 'IRIDIUM 161' number=43478 epoch=2019-02-23T06:25:35Z
[ 3062.67696436 -4537.7038142   4554.02354446]
Latitude: 39deg 59' 04.2"
Longitude: -12deg 29' 18.0"
Elevation (m): 751740
Traceback (most recent call last):
  File "skyfield-test.py", line 15, in <module>
    difference = sat_pos_geocentric - bluffton
  File "C:\Users\Kahn.non-free\.virtualenvs\inreach-kml-util-4s0grwMn\lib\site-packages\skyfield\positionlib.py", line 63, in __sub__
    p = self.position.au - body.position.au
AttributeError: 'Topos' object has no attribute 'position'

[brandon-rhodes]

The position of a lat/lon location changes every second as the Earth rotates. To produce a position from a Topos, try calling its .at(t) method with a particular time.

Where were you reading in the docs? If you let me know, I'll go to that spot and add a description clarifying this point!

[Kahn]

Thanks for coming back so quickly @brandon-rhodes the example I was following was at https://rhodesmill.org/skyfield/earth-satellites.html#generating-a-satellite-position but now you mention the .at(t) that makes perfect sense. I'll give it a go tonight

[Kahn]

@brandon-rhodes in the end I found out what I was doing wrong on a couple of accounts. Here is the working function I was working towards. Turns out when iterating on satellites you don't actually have an EarthSatellite.

def iridium_positions(utc_dt, lat_deg, lon_deg, elevation_m):

    stations_url = 'http://celestrak.com/NORAD/elements/iridium-NEXT.txt'
    satellites = load.tle(stations_url)

    for sat in satellites:
        if "IRIDIUM" in str(sat):
            elevation = str(elevation_m).split(" ")
            ts = load.timescale()
            t = ts.utc(utc_dt)
            geocentric = satellites[sat].at(t)
            # print("geocentric.position.km:", geocentric.position.km)

            subpoint = geocentric.subpoint()
            # print('Latitude:', subpoint.latitude)
            # print('Longitude:', subpoint.longitude)
            # print('Elevation (m):', int(subpoint.elevation.m))

            team = Topos(lat_deg, lon_deg, elevation_m=float(elevation[0]))
            difference = satellites[sat] - team
            # print('Difference:', difference)

            topocentric = difference.at(t)
            # print('Topocentric:', topocentric.position.km)

            alt, az, distance = topocentric.altaz()

            if alt.degrees > 9:
                print('The satellite {} is above the horizon').format(sat)
                print('Altitude Degrees:', alt.degrees)
                print('Azimuth Degrees:', az.degrees)
                print('Distance:', distance.km)`

[brandon-rhodes]

What do you mean by "you don't actually have an EarthSatellite"? Is sat a different sort of object?

In any case, I'm glad you figured everything out and got a working script! Thanks for posting the working code, it may help future readers who find this discussion through a web search.

[Kahn]

Yes, sat is just a string. The debugging that put me on the right track was printing the iterated output.

for sat in satellites:
    print(sat) # Gets a satellite name string like "IRIDIUM 123"
    print(satellites[sat]) # Gets an EarthSatellite by name