skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

Subpoints from Array of Geocentric Satellite Positions #230

Closed Tucktuckgoose closed 5 years ago

Tucktuckgoose commented 5 years ago

I am relatively new to both python and skyfield, and I am messing around to see what the capabilities are. I have successfully been able to calculate ISS coordinates in Lat, Long, Alt for a single point in time, but I am trying to write a program to generate a csv of those positions across a range of times. Here is my code up to this point.


from skyfield.api import Topos, load
import csv
from datetime import datetime

time_now = datetime.utcnow()

#Set the number of minutes on either side of time now that you want data for
time_buffer = 15
#Set the time buffer step in seconds
time_buffer_step = 5

#Read in the most up to date TLE for the ISS, and load it
stations_url = 'http://celestrak.com/NORAD/elements/stations.txt'
satellites = load.tle(stations_url)
satellite = satellites['ISS (ZARYA)']

with open('iss_location.csv', mode='w') as satellite_file:
    #Open a new CSV for writing
    satellite_write = csv.writer(satellite_file,delimiter=',',quotechar='"', quoting=csv.QUOTE_MINIMAL)

    #Create Headers
    satellite_write.writerow(['LAT','LONG','ALT'])

    #Calculate the buffer time
    time_low_seconds = time_now.second - (time_buffer * 60)
    time_high_seconds = time_now.second + (time_buffer * 60)

    #Set up the range of times to sample from
    ts =load.timescale()
    time_range = ts.utc(time_now.year, time_now.month, time_now.day, time_now.hour, time_now.minute, range(time_low_seconds, time_high_seconds, time_buffer_step))

    geocentric = satellite.at(time_range)

    for geo in geocentric:
        subpoint = geo.subpoint()
    #    satellite_write.writerow([subpoint.latitude.degrees,subpoint.longitude.degrees,int(subpoint.elevation.km)])

At the line subpoint = geo.subpoint() I get the error:

File "C:\Python27\lib\site-packages\skyfield\positionlib.py", line 465, in subpoint raise ValueError("you can only ask for the geographic subpoint" ValueError: you can only ask for the geographic subpoint of a position measured from Earth's center

Is there a problem with creating subpoints from an array, or am I doing something wrong?

Thanks and Merry Christmas!

brandon-rhodes commented 5 years ago

I'm not sure what happens if you try to loop over a geocentric object — you could add print(geo) to the loop so we can see what kind of object you got?

mattabern commented 5 years ago

I'm having basically the same error in a similar circumstance. When I run a print line in the loop it prints the object as <Geocentric position and velocity at date t>.

It looks like it's not handling the time array well. If I instead run with a single time instead of an array of times, it solves just fine.

The .position method works just fine with the time array.

mattabern commented 5 years ago

Well, I can force the issue with something like geo.center=399. Then it runs without issue.

What is the purpose of the if self.center != 399: check on line 464 of positionlib.py anyway? Why 399?

ghost commented 5 years ago

399 is the Earth's center (the geocenter) as defined by NAIF: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/naif_ids.html#Planets%20and%20Satellites

Update: my guess below is wrong-- line 84 of sgp4lib.py sets center = 399. So maybe try mapping over the list or, as Brandon suggested, printing to make sure you're actually getting back what you want.

THE PARAGRAPH BELOW IS WRONG:

I'm guessing the object satellite.at(time_range) returns is measured geocentrically, but geo.center isn't set to 399 for some reason. I think tweaking satellite.at would fix this.

brandon-rhodes commented 5 years ago

Happily, this appears to have been fixed during a recent refactoring! I can't reproduce it now with Skyfield master. Please re-open if this recurs in the future.