Skyfield should support Type 21 ephemeris segments

[JoshPaterson]

I have been using the Horizons telnet interface to generate binary SPK files for small solar system bodies. When I try to load them in skyfield I get this error:

File "skyfield\jpllib.py", line 179, in __new__
    .format(spk_segment.data_type))

ValueError: SPK data type 1 not yet supported segment

which seems to be referring to the data types discussed here: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html#Supported%20Data%20Types

Horizons is able to produce SPK files in a few different formats. The extended help for the file type selection prompt says this:

  [B]inary ... produce BINARY SPK file (.bsp). RECOMMENDED.
  [T]ext   ... produce ASCII transfer file (.xsp). DEPRECATED.
  1        ... produce legacy Type 1 binary format SPK file. DEPRECATED.
  -        ... move back to previous prompt
  ?        ... display a brief help file for this prompt
  ?!       ... display this extended help file
  X or Q   ... immediately disconnect from Horizons system

  * The initial default SPK format is "BINARY".  It is recommended this option
    be used.

...

  * The '1' response produces a Type 1 binary SPK file for use with older
    software linked against version of the SPICE Toolkit library prior to
    N0065. It is recommended users update to more recent SPICE reader libraries
    and not use this legacy option.

  * The alternate ASCII text form is deprecated and no longer necessary for
    users who link their software against N0052 and later version of the SPICE
    Toolkit reader library. However, if selected anyway, users must then
    convert the resulting ASCII format file to binary by running the SPICE
    utility program 'tobin', or 'spacit' on their local machines. Note that
    Horizons first produces the file in binary format then converts it to
    the ASCII transfer format if requested.

This seems to be implying that there are regular binary SPK files, and also legacy type 1 binary SPK files, but files created with either of these options all produce the same error I included above. I've tried this for a number of the largest small bodies (ceres, pallas, eris, etc.), and Horizons seems to only be capable of producing type 1 SPK files.

Based on a look through skyfield.jpllib.py it seems that skyfield currently supports data types 2 & 3.

The problem could be that I'm using Horizons incorrectly, but if it's not are there any plans to support data type 1 SPK files? That would allow skyfield to be used for asteroids, dwarf planets other than pluto, and comets.

Thanks!

P.S. Another small issue is that after failing to load a type 1 SPK file, I tried to rename it and I got a dialog saying "The action can't be completed because the file is open in Python." Is the file not being properly closed if that exception is raised?

[brandon-rhodes]

My guess is that their text is trying to distinguish not between "legacy Type 1 files" and "non-legacy non-Type-1 files", but between two kinds of Type 1 file: one that's legacy and deprecated, and another that's fully modern and supported. So I don't think you're misusing HORIZONS but just encountering an ambiguity in their documentation.

I would indeed like to add support for Type 1 ephemerides to Skyfield! I'm not sure when that will be, but at some point I'll go try reading the spke01_ Fortran function referenced in the SPK documentation in the Type 1 section:

https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html#Type 1: Modified Difference Arrays

Yes, that exception makes it sound like the file was not closed when you tried to rename it. I suppose I should at some point add a close() method to SpiceKernel so that users with a long-running Python script or an interactive session, and who are on operating systems that don't allow open files to be renamed, can easily close them. Let's keep this issue open until both of these features are added to Skyfield!

[JoshPaterson]

Thanks!

Also, I happened upon this script on the JPL site for automating access to Horizons and producing SPK files: ftp://ssd.jpl.nasa.gov/pub/ssd/SCRIPTS/smb_spk

In it's documentation, it has options for what file type to create, and one of them is this:

#            -2 create Type 21 binary SPK format. This new format will
#                become the default Nov 1, 2015 (what will be returned by 
#                the '-b' flag). In some situations, mostly planetary
#                encounters and over long periods of time, it can provide
#                a more accuracy, but requires linking against SPICE
#                Toolkit library N0065 or later, to access the new reader.

So it seems like the default binary SPK is type 21 (extended modified difference arrays), and the legacy binary SPK is type 1 (modified difference arrays). There is a fortran script spke21_ mentioned in the type 21 section here: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html#Type%2021:%20Extended%20Modified%20Difference%20Arrays

[brandon-rhodes]

Thanks for tracking down that bit of documentation! But one thing confuses me a bit: if the new binary format is type 21, then why did the error you received from Skyfield mention SPK data type 1? I had thought that the data type was read directly from the SPK file, and would say 21 if that was the value of the data type field in the file itself.

[JoshPaterson]

I set about double checking this by downloading new files from Horizons. The prompt for data type options is:

SPK text transfer format  [ YES, NO, ? ] :    

And the help says:

  [B]inary ... produce BINARY SPKfile (.bsp). RECOMMENDED.
  [T]ext   ... produce ASCII transfer file (.xsp). DEPRECATED.
  1        ... produce legacy Type 1 format binary file. DEPRECATED.
  -        ... move back to previous prompt
  ?        ... display this brief help file
  ?!       ... display an extended help file for this prompt
  X        ... immediately disconnect from Horizons system

Based on this, I expected that if I said 'yes' or 'text' I would get the ascii file, if I said 'no' or 'binary' I would get the type 21 binary file, and if I said '1' I would get the type 1 binary file.

But this is not what happens. Horizons gives me the exact same type 1 binary file whether I say 'no', 'binary', or '1'. This is why the error I got from Skyfield said both files were data type 1, they were actually the same file.

I then discovered there is an undocumented option '21' that produces a type 21 binary file. The error from Skyfield for this file does say type 21. I also tried numbers 2-20 to see if other data types were possible, but Horizons doesn't accept them.

So long story short, I think there's a bug in Horizons that makes the 'binary' option return the deprecated type 1 binary file, but the undocumented option '21' produces the type 21 binary file.

[brandon-rhodes]

Thank you for doing this additional investigation — it makes things much clearer! Are you going to report the problem to the HORIZONS folks to learn whether they're planning on making Type 21 the default, or at least documenting that the option exists?

[JoshPaterson]

I sent an email to Jon Giorgini at JPL, I believe he maintains Horizons. I'll let you know if I hear back.

I also found a relevant item on the Horizons news page:

Sep 23, 2015:
    -- Based on feedback from users, the switch to the new small-body SPK 
        default format will occur no sooner than November 1. Those using a 
        version of the 'smb_spk' automation script older than 2004-Aug-09 
        (version 1.3 or earlier) will need to update that also. The current 
        version 2.0 of the script is available here:

              ftp://ssd.jpl.nasa.gov/pub/ssd/SCRIPTS/smb_spk

        Remember to change the first line of the script to match the path 
        output by the UNIX/Linux command 'which expect' on your local system. 

        If you are using your own automation, you may need to adjust that once
        the change takes effect. The prompt string that currently reads:

          ' SPK text transfer format  [ YES, NO, ? ] : '

        ... will become:

          ' SPK file format  [ TRANSFER, BINARY, ? ] : '

        Note that the JPL script supplies an I/O model code to Horizons
        which allows the script to continue to work even if Horizons is
        changed in the future, so it may be more robust to use the provided 
        scripts.

[brandon-rhodes]

Thank you! So it looks like a transition they're taking slowly because scripts will break that expect the old format. That does make sense.

[JoshPaterson]

I heard back from Jon Giorgini, and he said that the users who needed it the old way have since fixed their scripts, and he changed it so option 'binary' now produces a type 21 file.

[brandon-rhodes]

Nice! So the focus of contributors to Skyfield should now be on type 21, not old type 1 files, correct?

[JoshPaterson]

Yes, that's correct.

[brandon-rhodes]

I'm marking this as contributor ready since I won't have time to do this soon, but someone else could do it by studying the existing code for bsp files and adding a new case for Type 21. The author would need to read the section "Type 21: Extended Modified Difference Arrays" of this document:

https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html#Type%2021:%20Extended%20Modified%20Difference%20Arrays

And they could probably start by porting the code fairly directly from the SPICE spke21_() function mentioned in that document.

[whiskie14142]

@brandon-rhodes I wrote a python module named 'spktype21' that reads binary SPK file of data_type 21 and computes positions and velocities of celestial small object.

You can get spktype21 from following links: https://github.com/whiskie14142/spktype21 https://pypi.org/project/spktype21/

Also, you can install it by pip command:

pip install spktype21

[whiskie14142]

@brandon-rhodes Please remind that function 'compute_type21' of spktype21 returns velocities in kilometers per second, in contrast, 'compute_and_differentiate' of jplephem returns velocities in kilometers per day.

[brandon-rhodes]

Thanks for pointing out your implementation, @whiskie14142! I've added this to my list of things to go read when I next sit down to work on Skyfield.

[davidmikolas]

@whiskie14142 I'd used your script earlier https://space.stackexchange.com/a/19263 and it worked nicely, thank you again! I see you've recently made some improvements.

[brandon-rhodes]

Skyfield’s documentation now officially recognizes @whiskie14142’s library as support for these ephemeris segments:

https://rhodesmill.org/skyfield/planets.html#third-party-libraries-for-other-ephemeris-formats

Hopefully someday there will be time to integrate the support into Skyfield’s core, but at least for now it’s fully documented.

[JoshPaterson]

I'm putting this link here so I can find it more easily in the future:

To generate a type 1 or type 21 SPK file for a single small body through a website GUI rather than the telnet interface: https://ssd.jpl.nasa.gov/x/spk.html

[stewi2]

Not quite sure what's going on, but whatever I do, Horizons will only give me an SPK file centered at the sun (10), rather than the solar system barycenter. Here's an example API request (same happens via the web form):

https://ssd.jpl.nasa.gov/api/horizons.api?MAKE_EPHEM=YES&COMMAND=%27DES%3D3727039%3B%27&EPHEM_TYPE=SPK&CENTER=%27500%400%27&OBJ_DATA=NO&START_TIME=2020-01-01&STOP_TIME=2030-12-31

And this is summary of the first segment (via print(kernel.segments[0])):

2459928.50..2459951.50  Sun (10) -> Unknown Target (3727039) data_type=21

Any idea what to do here? Right now the example code errors out with ValueError: Invalid Target and/or Center, since it's looking for a segment with center=0. When I change the call to look for center=10, it works, but that provides wrong results when I actually use it in Skyfield.

[brandon-rhodes]

@stewi2 — Searching the Skyfield code base, I can't find the phrase Invalid Target, so I'm not sure how to investigate your problem. Could you share a full traceback showing module name and line numbers? Thanks!

[stewi2]

It's coming from spktype21 when it's trying to find a segment that matches both the center and the object identifier. I can fix this by changing the 0 in this line to 10, but then things get messed up because self.center = 0 in the Type21Object class. And changing that to 10 makes it complain that I can't observe from a non-barycentric frame.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
[<ipython-input-136-fe62ff9b7925>](https://localhost:8080/#) in <module>
----> 1 alt, az, distance, alt_rate, az_rate, range_rate = raleigh.at(ts.utc(2022,12,18)).observe(rn35).apparent().frame_latlon_and_rates(coords)
      2 
      3 rate = AngleRate._from_radians_per_day(
      4     math.sqrt(alt_rate.radians.per_day**2 + az_rate.radians.per_day**2))
      5 
[/usr/local/lib/python3.8/dist-packages/skyfield/positionlib.py](https://localhost:8080/#) in observe(self, body)
    698 
    699         """
--> 700         p, v, t, light_time = body._observe_from_bcrs(self)
    701         astrometric = Astrometric(p, v, t, self.target, body.target)
    702         astrometric._ephemeris = self._ephemeris

[/usr/local/lib/python3.8/dist-packages/skyfield/vectorlib.py](https://localhost:8080/#) in _observe_from_bcrs(self, observer)
    104                              ' but this vector has the center {0}'
    105                              .format(self.center_name))
--> 106         return _correct_for_light_travel_time(observer, self)
    107 
    108     def geometry_of(self, other):

[/usr/local/lib/python3.8/dist-packages/skyfield/vectorlib.py](https://localhost:8080/#) in _correct_for_light_travel_time(observer, target)
    241     cvelocity = observer.velocity.au_per_d
    242 
--> 243     tposition, tvelocity, gcrs_position, message = target._at(t)
    244 
    245     distance = length_of(tposition - cposition)

[<ipython-input-135-013e0e6b16a4>](https://localhost:8080/#) in _at(self, t)
     11     def _at(self, t):
     12         k = self.kernel
---> 13         r, v = k.compute_type21(0, self.target, t.whole, t.tdb_fraction)
     14         return r / AU_KM, v / AU_KM, None, None
     15 

[/usr/local/lib/python3.8/dist-packages/spktype21.py](https://localhost:8080/#) in compute_type21(self, center, target, jd1, jd2)
    131                 return result[0:3], result[3:]
    132 
--> 133         self.mda_record, self.mda_lb, self.mda_ub = self.get_MDA_record(eval_sec, target, center)
    134         self.mda_record_exists = True
    135 

[/usr/local/lib/python3.8/dist-packages/spktype21.py](https://localhost:8080/#) in get_MDA_record(self, eval_sec, target, center)
    168                 matched.append(segment)
    169         if len(matched) == 0:
--> 170             raise ValueError('Invalid Target and/or Center')
    171         if eval_sec < matched[0].start_second or eval_sec >= matched[-1].end_second:
    172             raise ValueError('Invalid Time to evaluate')

ValueError: Invalid Target and/or Center

[brandon-rhodes]

Alas, as I don't know the internals of how the spktype21.py file works, I'm not sure I can offer advice. In Skyfield, a segment that starts at the Sun can be used if you do something like:

sun = eph['Sun']
obj = sun + segment

— where eph is something like DE421 that gives you the position of the Sun over time, and segment is the position that leads from the Sun out to the object of interest. But I'm not sure how to apply that to your code without an example script that illustrates what you're dealing with, alas.

[stewi2]

Yes - that worked - thank you! I still have to get used to the coordinate transformation logic :) It might be worth updating the example here: https://rhodesmill.org/skyfield/planets.html#type-1-and-type-21-ephemeris-formats, since it assumes that Horizons returns SPK files centered at 0, when they're really centered at 10. But as you point out I can also easily transform the resulting vector. I'm wondering if something changed in Horizons.