skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.4k stars 209 forks source link

Creating a Skyfield Body from JPL Horizons data for Spacecraft (not comets/asteroids)? #350

Closed zleffke closed 4 years ago

zleffke commented 4 years ago

Hello! This isn't really an issue, just a request for guidance, my apologies if this should be posted elsewhere. I would like to learn how to create a skyfield 'body' from spacecraft ephemeris (time series of state vectors) downloaded from JPL HORIZONS.

I can successfully query the JPL HORIZONS database to obtain a time series of state vector data for a spacecraft (LRO for this example) with the Solar System Barycenter as the reference origin. What I would like to be able to do is replicate the functionality for computing observations to the target spacecraft, just like one would do for computing pointing angles to something loaded from the de421.bsp file (for example). I have been able to create individual Skyfield Barycentric positions from the HORIZONS data, but what I would like to do is create an actual LRO 'body' so that the rest of the code would follow the same process of using the 'observe' method. In addition to pointing antennas, I would also like to do things like figure out when LRO is occulted by the Moon and thus compute access windows, so just computing pointing angles to the moon is not sufficient.

I've looked into various methods for downloading or generating binary SPKs from various online references (including some of the issues on this repo) but the SPKs are either of the wrong type for skyfield (not type 2 or 3), mkspk from NAIF doesn't produce type 2 or 3, and the links for autogenerating SPKs do not work for spacecraft (asteroids and comets only). I suspect that generating an SPK of type 2 or 3, using SpicePy and Chebyshev regression based on time series state vector data from HORIZONS might work, and then could be loaded by Skyfield, but that feels a bit round about (and I'm not sure I actually know how to do that....maybe something with numpy/scipy for the regression part to get the coefficients?)

Ultimately I would like to be able to do something like:

gs.at(t).observe(lro)

where the LRO body is based on the HORIZONS state vector data.

Any guidance on how to make spacecraft 'bodies' in Skyfield from HORIZONS data? is this something that could perhaps be accomplished using jplephem?

Thanks in advance! Skyfield rocks!

brandon-rhodes commented 4 years ago

Skyfield rocks!

Thank you — I'm glad that you are enjoying Skyfield so far!

What type number of SPK are you able to acquire for spacecraft, and from which service? The mostly likely answer is that we should add support for that SPK type to jplephem and then explain to Skyfield that it can pass a date in and get a position back out, like it does for planets with their Type 2 and Type 3 files.

zleffke commented 4 years ago

Hi Brandon, I have yet to directly find an existing SPK for LRO. That said it appears there are a number of potential sources for this sort of thing (not just LRO), and I'm sure you are aware of most (if not all) of these. I'll attempt to itemize what I've explored here in case its useful for others.

  1. https://naif.jpl.nasa.gov/pub/naif/. This seems to have a number of various spacecraft listed. It is not 100% comprehensive and identical to HORIZONS. LRO SPKs for example are not available at this link, but you can get current LRO data from HORIZONS. I have not extensively examined the types available here, but suspect it is a mix. Other references for this site show LRO as 'partially' supported (CK/FK files are there). Juno BSP juno_pred_orbit.bsp yields type 1. Issue 157 discusses Type 1 support (though issue is specifically marked for Type 21 contributor ready support). STEREO-A_merged.bsp yields type 13. mro_psp.bsp yields type 1. So it looks like a mix.
  2. https://ssd.jpl.nasa.gov/x/spk.html. Web based utility for autogenerating SPKs for 'small body' objects such as asteroids and comets from HORIZONS. Does not generate SPKs for spacecraft (will throw an error if querying it for a spacecraft ID). Following the example in Issue 196 for Asteroid 2007 RH283 yields an SPK of type 21. Issue 157 already exists for adding type 21 support.
  3. Horizons Web Interface. Generates state vector data for spacecraft, but it does not directly generate binary SPKs. Only way to download data is in plain text form and is a time series of state vector data.
  4. Horizons Telnet Interface. Tried this as well, but like number 2 above it does not appear to be useful for generating SPKs. Pretty sure this is because spacecraft are not asteroids/comets where Binary SPK generation (apparently type 21) is possible.
  5. mkspk Utility from NAIF. From their description: "MKSPK is a NAIF Toolkit utility program that generates a new, or appends data to an existing, spacecraft or target body's ephemeris file in SPICE SPK format using one of the following SPK types: 5, 8, 9, 10, 12, 13,15, 17." I suspect there may be a way to use this along with data queried according to number 3 or 4 above to generate an SPK. Based on the description from the NAIF Spice Tutorials, tutorial 42_making_an_spk, slide 16, it would appear types 8, 9, 12, & 13 might be appropriate. This method however seems like a hack because the source data from JPL HORIZONS is probably already based on regression.....so you would be querying for values based on regression..to create regression coefficients...to create an SPK based on regression...to query for values....this is the method that seems roundabout (and perhaps introduces errors or inaccuracies with the 'double regression').
  6. SpiceyPy. Similar to number 5 but a pythonic method for using SPICE utilities. I was looking at this as an alternative to number 5 to generate SPKs of type 2 or 3 that would be accepted by Skyfield. Again though, seems like a hack to pull HORIZONS data (probably from regression data), generate a regression, then generate the SPK. roundabout again.

That about summarizes all the methods I've explored so far. Some of the BSPs for spacecraft from NAIF are in type 1 an 13 format (non-comprehensive spot checks, there could be other types), but NAIF doesn't have all the spacecraft that HORIZONS has. So doesn't look like there is anything new you haven't already referenced in other issues, unless the types mentioned in number 5 above are of interest (they seem to be relevant to spacecraft trajectories).

There doesn't seem to be a method to directly generate SPKs for spacecraft from HORIZONS. Is there a method for getting the source data from HORIZONS that they use to generate the state vectors for spacecraft? It might be mission dependent, but perhaps they already have the Chebysev coefficients (or similar) that could be converted to SPK/BSP format (by the way, whats the difference between SPK and BSP file extensions?).

Hopefully the list above is helpful. Not sure of a path forward for this unless you or others have any ideas. Ultimately, I think it will be really useful for folks to be able to track spacecraft not trackable by TLEs (not in Earth Orbit), so would like to keep pulling this thread. One example is future Amateur Radio missions that might involve the Lunar Gateway in the NRHO they are planning. Again, simply pointing at the Moon may not be sufficient for the NRHO, especially if using high gain microwave antennas with narrow beams. This of course assumes SPKs and/or Gateway ephemeris shows up in either NAIF or HORIZONS databases.

Thanks again for your response. Your general responsiveness and support is amazing.

-Zach, KJ4QLP

zleffke commented 4 years ago

I should also mention I've used astroquery to directly access HORIZONS from Python scripts based on this documentation. Useful (skips the file download/python import and parse steps), but ultimately the same as 3 and 4 in the previous comment.

brandon-rhodes commented 4 years ago

To (finally) follow up here: thanks for such a complete answer, listing all the approaches you've tried.

No, I don't know where HORIZONS gets its spacecraft data. I wish that its source code were available as well as the data behind it, but they seem to keep that all under wraps.

Right now Type 1 and Type 21 are supported through a pair of packages independent of Skyfield:

https://github.com/whiskie14142/spktype01

https://github.com/whiskie14142/spktype21

I suppose I should sit down sometime soon and see about getting that code imported so that it becomes a native feature.

brandon-rhodes commented 4 years ago

Using Type 1 and Type 21 ephemeris formats in Skyfield was recently documented on its "Planets" page, down near the bottom:

https://rhodesmill.org/skyfield/planets.html

While it requires third-party packages, it does seem to work, so I am going to close this issue for the moment. But I will add a long-term item in the TODO.rst file reminding contributors that someday we would like to vectorize those computations and bring them inline into Skyfield!