brandon-rhodes
Many time intervals for segments with same values
Hello
There are some BSP files I work that have many segments with same values, each with different time intervals.
For example:
<-->**k=SPK.open("chariklo_jpl20.bsp")**
<-->**print(k)**
....
2487340.50..2487468.50 Sun (10) -> Unknown Target (2010199)
2487468.50..2487596.50 Sun (10) -> Unknown Target (2010199)
2487596.50..2487724.50 Sun (10) -> Unknown Target (2010199)
2487724.50..2487852.50 Sun (10) -> Unknown Target (2010199)
2487852.50..2487980.50 Sun (10) -> Unknown Target (2010199)
2487980.50..2488108.50 Sun (10) -> Unknown Target (2010199)
2488108.50..2488236.50 Sun (10) -> Unknown Target (2010199)
2488236.50..2488364.50 Sun (10) -> Unknown Target (2010199)
2488364.50..2488428.50 Sun (10) -> Unknown Target (2010199)
2488402.50..2488434.50 Sun (10) -> Unknown Target (2010199)
<---> **k[10,2010199].compute(2457023.5)**
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-46-333d3990b7b5> in <module>()
----> 1 k[10,2010199].compute(2457023.5)
/home/altairjunior/anaconda/lib/python2.7/site-packages/jplephem/spk.pyc in compute(self, tdb, tdb2)
104 def compute(self, tdb, tdb2=0.0):
105 """Compute the component values for the time `tdb` plus `tdb2`."""
--> 106 for position in self.generate(tdb, tdb2):
107 return position
108
/home/altairjunior/anaconda/lib/python2.7/site-packages/jplephem/spk.pyc in generate(self, tdb, tdb2)
164 final_epoch = initial_epoch + interval_length * n
165 raise ValueError('segment only covers dates %.1f through %.1f'
--> 166 % (initial_epoch, final_epoch))
167
168 omegas = (index == n)
ValueError: segment only covers dates 2488402.5 through 2488434.5
So, it only gets the last segment and does not compute for instant in the middle
Is there any way to avoid this?
Oh, wow, neat — you get to use more complicated *.bsp files than the ones I have seen so far!
Yes, there is a way around. The k[…] lookup convention is just a convenience: a shorthand for saying k.segments[i] where i is the index into the list of segments. To grab a specific segment you could do something like:
k.segments[5].compute(2457023.5)
To select the correct segment automatically (though this won’t be terribly efficient for, say, millions of data points if you do this with every one):
def segment_for(k, center, target, jd):
for segment in k.segments:
if segment.center == center and segment.target == target and segment.start_jd <= jd <= segment.end_jd:
return segment
raise ValueError('no segment matches')
It would be more efficient, of course, to pre-sort the list of segments by center-and-target so that you could grab the correct segment list and then only have to compare the dates as you ran through it, since you would know that the center and target were already correct. And there is probably a faster way than iteration to correctly choose the segment for a given date, especially if the segments all have the same length because then you could just divide to produce an integer index.
But for now, try something like the above to quickly get yourself up and running, and let me know whether it works! Maybe a future version of jplephem should include a more sophisticated way to choose segments when there are several for the same objects.
Do you know why this ephemeris doesn’t just have a single big segment, but instead has several little ones? I could not find a copy of it by searching for the filename on Google.
Hello
Thank you, it really helped
First, I want to point in your code that instead of "k.center" it is "segment.center". The same goes for "k.target"
I don't know the sources of these bsp files since it's a collaborator who maintains our server. But I know that in this link you can download the "ph15.bsp" file that have the same issue. http://josselin.desmars.free.fr/ph15/ It is an ephemeris for the Saturnian satellite Phoebe made by a french collaborator.
I don't know why they are constructed like this, since I don't know how they are constructed. I though it was usual.
Thanks, I fixed the code in my comment!
Let’s keep this issue open. When I get the chance someday, I will look into including in jplephem an algorithm for efficiently accessing ephemerides that have more than one segment for the same center and target — and thanks for the link, I will use ph15.bsp as my example!
Here is a JPL file with two segments per body:
ftp://ssd.jpl.nasa.gov/pub/eph/satellites/bsp/jup357.bsp
SPICE kernel file 'jup357.bsp' has 18 segments
JD 2450464.50 - JD 2524606.50 (1997-01-15 through 2200-01-13)
5 -> 599 JUPITER BARYCENTER -> JUPITER
5 -> 516 JUPITER BARYCENTER -> METIS
5 -> 515 JUPITER BARYCENTER -> ADRASTEA
5 -> 514 JUPITER BARYCENTER -> THEBE
5 -> 505 JUPITER BARYCENTER -> AMALTHEA
5 -> 504 JUPITER BARYCENTER -> CALLISTO
5 -> 503 JUPITER BARYCENTER -> GANYMEDE
5 -> 502 JUPITER BARYCENTER -> EUROPA
5 -> 501 JUPITER BARYCENTER -> IO
JD 2378482.50 - JD 2450464.50 (1799-12-17 through 1997-01-15)
5 -> 599 JUPITER BARYCENTER -> JUPITER
5 -> 516 JUPITER BARYCENTER -> METIS
5 -> 515 JUPITER BARYCENTER -> ADRASTEA
5 -> 514 JUPITER BARYCENTER -> THEBE
5 -> 505 JUPITER BARYCENTER -> AMALTHEA
5 -> 504 JUPITER BARYCENTER -> CALLISTO
5 -> 503 JUPITER BARYCENTER -> GANYMEDE
5 -> 502 JUPITER BARYCENTER -> EUROPA
5 -> 501 JUPITER BARYCENTER -> IO
@jackschmidt — Interesting! I wonder if the before-1997 and after-1997 segments have different resolutions?
I think the segments have the same resolution, but different time periods. The moons have a different resolution than the planet. I've just started using your code and JPL's data today, so please forgive any mistakes.
>>> start,interval,coeffs=skyfield.api.load("jup357.bsp").segments[0].spk_segment.load_array()
>>> (start,interval,coeffs.shape)
(2450464.5, 1.0, (3, 74142, 11))
>>> start,interval,coeffs=skyfield.api.load("jup357.bsp").segments[9].spk_segment.load_array()
>>> (start,interval,coeffs.shape)
(2378482.5, 1.0, (3, 71982, 11))
@jackschmidt — With identical intervals, I’m puzzled why they didn’t just make one segment for each moon instead of two for each moon.
Skyfield currently has not automatic way to switch between segments when asked about different time periods. Could you use manual indexing into the segments list for now, to grab the moons you want to use? Let’s figure out what’s missing from the current Skyfield documentation to guide you in solving this problem, and I can go write up the answers you need!
I agree it is weird and honestly my next step may be to write a bsp file editor that pulls out the date range needed (possibly from multiple files) and creates a new (for me, smaller) file.
For Skyfield: I haven't found a general solution, but I used your code above, updated to the current naming conventions and argument types. I am mostly interested in times near the present, so my solution was to specify ts.now() as the target time for finding the segments.
There is a related annoyance: if I want to work with Europa, I need to load both a "de" and a "jup" bsp file, and I have to somewhat manually lookup the correct chain of vectors. The relationship is specifically: "to load a solar system barycenter to solar system body vector over a time range may require multiple bsp files". Perhaps there could be something similar to the class SPK that can handle multiple files simultaneously. I think solving this weird multiple segment issue is about as complicated as solving the multiple file version.
My code looks like this currently:
def segment_for(self, center, target, when):
for bsp in self.bsps:
for segment in bsp.segments:
(start,stop) = segment.time_range(self.ts)
if segment.center == center and segment.target == target and start.tt <= when.tt <= stop.tt:
return segment
raise ValueError('no segment matches')
def body(self, bodyName):
if bodyName in self.planets:
return self.planets[bodyName]
elif bodyName in self.jupiter:
return self.planets['jupiter barycenter'] + self.segment_for(
center = self.jupiter.decode('jupiter barycenter'),
target = self.jupiter.decode(bodyName),
when = self.ts.now() )
else:
return self.planets[bodyName + ' barycenter']
where awkwardly I have specific bsps in self.planets and self.jupiter, and then the list of all available bsp in self.bsps.
When trying to figure out the general solution, I thought about people who are interested in historical data too. Some of the files have a "_1600" version where even more time is split out. I could not muster the energy to figure out how to load earth.at(t) when the times t stretched across multiple segments or files.
@jackschmidt — In an open issue from a couple of years ago, we worked out how to extend one ephemeris with segments from another:
https://github.com/skyfielders/python-skyfield/issues/145#issuecomment-357561455
I wonder if you could load both ephemerides, then do something like:
self.planets.segments.extend(self.jupiter.segments[:9]) # note: only taking the first 9 segments
Then maybe all the normal Skyfield logic for connecting segments back to the barycenter would take over, and your awkward code could disappear?
I suppose I should either create an official super-ephemeris object which can combine ephemeris files, or else at least add the above maneuver to the docs so that folks don't have to run across that earlier issue to learn the trick.
my solution was to specify ts.now() as the target time for finding the segments.
That's a good temporary solution! And is more robust than my take-the-first-nine approach in the code snippet above. :slightly_smiling_face:
I could not muster the energy to figure out how to load earth.at(t) when the times t stretched across multiple segments or files.
And I have not, alas, mustered the energy either. A solution will have to do things like be able to take an array of times for which positions are desired, and send different times to different segments depending on the dates they cover.
With one small addition this works very well:
self.planets.codes = self.planets.codes.union(self.jupiter.codes)
This way the very convenient self.planets['europa'] works. I haven't extensively tested, but I suspect this will just work.
This handles the multiple files issue quite well. The multiple segments per body problem remains unsolved, but it does not affect me. :-)
In case it helps future readers, my setup code is:
self.planets = load('de421.bsp')
self.jupiter = load('ftp://ssd.jpl.nasa.gov/pub/eph/satellites/bsp/jup357.bsp')
self.ts = load.timescale()
now = self.ts.now()
self.planets.segments.extend( s for s in self.jupiter.segments if s.time_range(self.ts)[1].tt > now.tt )
self.planets.codes = self.planets.codes.union(self.jupiter.codes)
and now self.planets works like planets and eph in the documentation, but supports a few moons of jupiter as well the other major solar system bodies.
Thanks for trying that idea out! I'll see if I can get the docs updated soon, to handle the case where a user can survive on a single time period rather than several overlapping segment time periods.
Another JPL file with multiple segments per body is de441.bsp on
Note that DE441 is split in two files on
- <https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de441_part-1.bsp
- <https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de441_part-2.bsp
$ python -m jplephem spk de441.bsp
File type DAF/SPK and format LTL-IEEE with 28 segments:
-3100015.50..2440432.50 Type 2 Venus Barycenter (2) -> Venus (299)
-3100015.50..2440432.50 Type 2 Mercury Barycenter (1) -> Mercury (199)
-3100015.50..2440432.50 Type 2 Earth Barycenter (3) -> Earth (399)
-3100015.50..2440432.50 Type 2 Earth Barycenter (3) -> Moon (301)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Sun (10)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Pluto Barycenter (9)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Neptune Barycenter (8)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Uranus Barycenter (7)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Saturn Barycenter (6)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Jupiter Barycenter (5)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Mars Barycenter (4)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Earth Barycenter (3)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Venus Barycenter (2)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Mercury Barycenter (1)
2440432.50..8000016.50 Type 2 Venus Barycenter (2) -> Venus (299)
2440432.50..8000016.50 Type 2 Mercury Barycenter (1) -> Mercury (199)
2440432.50..8000016.50 Type 2 Earth Barycenter (3) -> Earth (399)
2440432.50..8000016.50 Type 2 Earth Barycenter (3) -> Moon (301)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Sun (10)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Pluto Barycenter (9)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Neptune Barycenter (8)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Uranus Barycenter (7)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Saturn Barycenter (6)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Jupiter Barycenter (5)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Mars Barycenter (4)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Earth Barycenter (3)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Venus Barycenter (2)
2440432.50..8000016.50 Type 2 Solar System Barycenter (0) -> Mercury Barycenter (1)
$ python -m jplephem spk de441_part-1.bsp
File type DAF/SPK and format LTL-IEEE with 14 segments:
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Mercury Barycenter (1)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Venus Barycenter (2)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Earth Barycenter (3)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Mars Barycenter (4)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Jupiter Barycenter (5)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Saturn Barycenter (6)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Uranus Barycenter (7)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Neptune Barycenter (8)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Pluto Barycenter (9)
-3100015.50..2440432.50 Type 2 Solar System Barycenter (0) -> Sun (10)
-3100015.50..2440432.50 Type 2 Earth Barycenter (3) -> Moon (301)
-3100015.50..2440432.50 Type 2 Earth Barycenter (3) -> Earth (399)
-3100015.50..2440432.50 Type 2 Mercury Barycenter (1) -> Mercury (199)
-3100015.50..2440432.50 Type 2 Venus Barycenter (2) -> Venus (299)
$ python -m jplephem spk de441_part-2.bsp
File type DAF/SPK and format LTL-IEEE with 14 segments:
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Mercury Barycenter (1)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Venus Barycenter (2)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Earth Barycenter (3)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Mars Barycenter (4)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Jupiter Barycenter (5)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Saturn Barycenter (6)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Uranus Barycenter (7)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Neptune Barycenter (8)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Pluto Barycenter (9)
2440400.50..8000016.50 Type 2 Solar System Barycenter (0) -> Sun (10)
2440400.50..8000016.50 Type 2 Earth Barycenter (3) -> Moon (301)
2440400.50..8000016.50 Type 2 Earth Barycenter (3) -> Earth (399)
2440400.50..8000016.50 Type 2 Mercury Barycenter (1) -> Mercury (199)
2440400.50..8000016.50 Type 2 Venus Barycenter (2) -> Venus (299)
