Using skyfield to get small bodies such as Chiron

[emma88888888]

Hi there,

I remember seeing that you were planning to improve access to small bodies, such as Chiron, posted a few years back. Have you completed it yet?

Really keen to use this

[brandon-rhodes]

I started working on comets and minor planets again last week, and am planning to wrap it up this coming week! Do you know where you prefer to pull its data from — whether from the Minor Planet Center, or another data source?

[emma88888888]

I was hoping to use the Horizon data. I get a "ValueError: SPK data type 21 not yet supported segment" when trying to load the file

On Mon, Jun 8, 2020 at 12:04 PM Brandon Rhodes notifications@github.com wrote:

I started working on comets and minor planets again last week, and am planning to wrap it up this coming week! Do you know where you prefer to pull its data from — whether from the Minor Planet Center, or another data source?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/387#issuecomment-640299991, or unsubscribe https://github.com/notifications/unsubscribe-auth/APU2VXMP4GSWATPIUKXGHZLRVQTHZANCNFSM4NXXSQOQ .

[brandon-rhodes]

I was hoping to use the Horizon data. I get a "ValueError: SPK data type 21 not yet supported segment" when trying to load the file …

Ah, so in your case you are interested in a service that has already produced an SPK file.

There is a third-party package that supports Type 21 with Skyfield, but I have never had the chance to try it out myself and integrate it into Skyfield itself — if you try it out, please let me know about your experience! Here’s the link:

https://github.com/whiskie14142/spktype21

[emma88888888]

I tried it and failed unfortunately. See: https://github.com/whiskie14142/spktype21/issues/1

On Mon, Jun 8, 2020 at 12:26 PM Brandon Rhodes notifications@github.com wrote:

I was hoping to use the Horizon data. I get a "ValueError: SPK data type 21 not yet supported segment" when trying to load the file …

Ah, so in your case you are interested in a service that has already produced an SPK file.

There is a third-party package that supports Type 21 with Skyfield, but I have never had the chance to try it out myself and integrate it into Skyfield itself — if you try it out, please let me know about your experience! Here’s the link:

https://github.com/whiskie14142/spktype21

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/387#issuecomment-640303292, or unsubscribe https://github.com/notifications/unsubscribe-auth/APU2VXJA3EMSNXIYBVIZGOLRVQV27ANCNFSM4NXXSQOQ .

[brandon-rhodes]

“commented 19 days ago” — Oh, I see you have already done some research on other approaches!

I’ll see if I can schedule some time tomorrow to check out the library myself, look at the error message you are getting, and see what might be required to get the feature adapted for Skyfield itself.

[emma88888888]

thanks!

On Mon, Jun 8, 2020 at 2:02 PM Brandon Rhodes notifications@github.com wrote:

“commented 19 days ago” — Oh, I see you have already done some research on other approaches!

I’ll see if I can schedule some time tomorrow to check out the library myself, look at the error message you are getting, and see what might be required to get the feature adapted for Skyfield itself.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/387#issuecomment-640321484, or unsubscribe https://github.com/notifications/unsubscribe-auth/APU2VXP3US3M46IQKY5GOALRVRBD5ANCNFSM4NXXSQOQ .

[brandon-rhodes]

@emma88888888 — I was successfully able to induce HORIZONS to produce a Type 21 ephemeris for Chiron for me (I had not used telnet in years, so retro!), and was able to get the same error you did. So I looked at which bodies the ephemeris supports:

$ python -m jplephem spk ci/wld23593.15

File type DAF/SPK and format LTL-IEEE with 1 segments:
INITIAL JD..FINAL JD   TYPE   CENTER -> TARGET
2459008.50..2459015.50  Type 21  Solar System Barycenter (0) -> Unknown Target (2002060)

So there are two problems with your script, it would seem:

1. You are asking about a position based on body 399, the Earth, but HORIZONS might be giving you an ephemeris without any segments that provide the position of Earth?
2. Instead of using the raw minor planet number, it looks like HORIZONS adds 2 million to such numbers to avoid having them collide with other identification numbers, so it calls Chiron by the number 2002060 instead.

When using the numbers displayed above, I ask for:

position, velocity = chironData.compute_type21(0, 2002060, 2459009.3008)

— I get a position and velocity.

If this indeed solves the mystery of the error you have been seeing on your end (let me know how the ephemeris it gave me compares with yours), I will look at how big a project it might be to get this SPKType21 library working with the rest of Skyfield.

[brandon-rhodes]

@emma88888888 — It looks like it would take a few days of work to incorporate the spktype21 module into the underlying jplephem library and then teach Skyfield how to use it, and so I am not likely to be able to add official support within the next couple of weeks. But I have at least worked out a quick way to use the spktype21 third-party library as it currently stands with Skyfield:

from skyfield.api import load
from skyfield.constants import AU_KM
from skyfield.vectorlib import VectorFunction
from spktype21 import SPKType21

class Type21Object(VectorFunction):
    def __init__(self, kernel, target):
        self.kernel = kernel
        self.center = 0
        self.target = target

    def _at(self, t):
        k = self.kernel
        r, v = k.compute_type21(0, self.target, t.whole, t.tdb_fraction)
        return r / AU_KM, v / AU_KM, None, None

ts = load.timescale(builtin=True)
t = ts.utc(2020, 6, 9)

eph = load('de421.bsp')
earth = eph['earth']

kernel = SPKType21.open('ci/wld23593.15')
chiron = Type21Object(kernel, 2002060)

ra, dec, distance = earth.at(t).observe(chiron).radec()
print(ra)
print(dec)

This code is overly simple — it only works with objects at the Solar System barycenter — but it displays coordinates for Chiron that agree with those that HORIZONS produces for me:

00h 27m 38.99s
+05deg 57' 08.9"

If you could try this code out and confirm, I could at least add the above class to Skyfield’s documentation for other folks who need to use Type 21 ephemerides. Thanks for any feedback you can provide!

[emma88888888]

It worked perfectly. I just had to update to the latest skyfield for the improved time component.

thank you!

On Tue, Jun 9, 2020 at 8:03 AM Brandon Rhodes notifications@github.com wrote:

@emma88888888 https://github.com/emma88888888 — It looks like it would take a few days of work to incorporate the spktype21 module into the underlying jplephem library and then teach Skyfield how to use it, and so I am not likely to be able to add official support within the next couple of weeks. But I have at least worked out a quick way to use the spktype21 third-party library as it currently stands with Skyfield:

from skyfield.api import load

from skyfield.constants import AU_KM

from skyfield.vectorlib import VectorFunction

from spktype21 import SPKType21

class Type21Object(VectorFunction):

def __init__(self, kernel, target):

    self.kernel = kernel

    self.center = 0

    self.target = target

def _at(self, t):

    k = self.kernel

    r, v = k.compute_type21(0, self.target, t.whole, t.tdb_fraction)

    return r / AU_KM, v / AU_KM, None, None

ts = load.timescale(builtin=True)

t = ts.utc(2020, 6, 9)

eph = load('de421.bsp')

earth = eph['earth']

kernel = SPKType21.open('ci/wld23593.15')

chiron = Type21Object(kernel, 2002060)

ra, dec, distance = earth.at(t).observe(chiron).radec()

print(ra)

print(dec)

This code is overly simple — it only works with objects at the Solar System barycenter — but it displays coordinates for Chiron that agree with those that HORIZONS produces for me:

00h 27m 38.99s

+05deg 57' 08.9"

If you could try this code out and confirm, I could at least add the above class to Skyfield’s documentation for other folks who need to use Type 21 ephemerides. Thanks for any feedback you can provide!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/387#issuecomment-640856599, or unsubscribe https://github.com/notifications/unsubscribe-auth/APU2VXJERIK3S6FNUROC2VLRVU7ZXANCNFSM4NXXSQOQ .

[brandon-rhodes]

@emma88888888 — I'm glad it worked! I'll hopefully have time to get the documentation updated tomorrow to share this shortcut more widely.

[xmichaelx]

This is awesome development, thank you.

I'm testing this on 541 Deborah. I've generated spice kernel from Asteroid & Comet SPK File Generation Request site with following settings:

For comparison I'm using SkyBot service search (epoch embedded in url).

The code is copied from above example:

from skyfield.api import load
from skyfield.constants import AU_KM
from skyfield.vectorlib import VectorFunction
from spktype21 import SPKType21

class Type21Object(VectorFunction):
    def __init__(self, kernel, target):
        self.kernel = kernel
        self.center = 0
        self.target = target

    def _at(self, t):
        k = self.kernel
        r, v = k.compute_type21(0, self.target, t.whole, t.tdb_fraction)
        return r / AU_KM, v / AU_KM, None, None

ts = load.timescale(builtin=True)
t = ts.utc(2020,1,15,16,43,30) # epoch should be the same as in SkyBot search
print(t)
eph = load('de421.bsp')
earth = eph['earth']
kernel = SPKType21.open('2000541.bsp')
deborah = Type21Object(kernel, 2000541)
ra,dec,dist=earth.at(t).observe(deborah).radec()
print(ra,dec)

I see only small difference:

skyfield reports 07h 38m 11.10s +17deg 28' 22.5" skybot reports 07h 38m 11.01s +17deg 28' 22.616"

Is this something expected and can be interpreted as validation for this example? I'm adding script and exported kernel as attachment for reproducibility. deborah_test.zip

[brandon-rhodes]

@xmichaelx — Yes, I think that validates the example. A small discrepancy of a few arcseconds can be difficult to diagnose without access to the source code of both tools. Otherwise you have to go through every possible difference (Solar system ephemeris used, for example) and even then might never hit on the right tweak that makes them finally agree.

[brandon-rhodes]

I have added documentation to Skyfield that outlines the solution I shared above. Hopefully this will help folks until such time as the support can be rewritten as part of Skyfield. The documentation should appear on the web site when I release the next version of Skyfield later this week!

[emma88888888]

Hi Brandon,

I've just been working on an app that uses skyfield and thought you might really like it. Here's the link:

https://play.google.com/store/apps/details?id=com.blueskywoodpigeon.astromap

Cheers,

Femke

On Wed, Jun 10, 2020 at 2:17 PM Brandon Rhodes notifications@github.com wrote:

@emma88888888 https://github.com/emma88888888 — I'm glad it worked! I'll hopefully have time to get the documentation updated tomorrow to share this shortcut more widely.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/387#issuecomment-641678735, or unsubscribe https://github.com/notifications/unsubscribe-auth/APU2VXL5GJ5YSL4YTUDHYADRV3ULPANCNFSM4NXXSQOQ .

[brandon-rhodes]

Interesting, I do not remember ever seeing the idea of projecting constellation outlines on a map to make it visually clear where the stars stand above the Earth! Thanks for sharing the link, and I am glad you were able to make progress on your project using Python and Skyfield.