Vectorization of methods for minor bodies propagation

[xmichaelx]

This is basically copy-paste attempt of what @JoshPaterson did in #511

Following code shows it (MPCORB.DAT has header removed manually):

from skyfield.api import load
from skyfield.data import mpc
from skyfield.data.mpc import _mpcorb_orbits
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from numpy import arange, linspace, rollaxis
import matplotlib.pyplot as plt

ts = load.timescale()

with load.open('MPCORB.DAT') as f:
    rows = mpc.load_mpcorb_dataframe(f)
    bad_orbits = rows.semimajor_axis_au.isnull()
    rows = rows[~bad_orbits]
    rows = rows.set_index('designation', drop=False)
    comets = _mpcorb_orbits(rows.iloc[0:20], ts, GM_SUN)

    t = ts.utc(linspace(2000, 2030, 100))
    position = comets.at(t).position.au

    plt.rcParams['legend.fontsize'] = 10

    fig = plt.figure()
    ax = fig.gca(projection='3d')

    for i, (x, y, z) in enumerate(rollaxis(position, 1)):
        ax.plot(x, y, z, label=rows.designation[i])

    ax.legend()

    plt.show()

Produces following plot:

This still needs work but I'd like to use PR for discussion and to show some nice tricks for speeding up epoch_packed parsing in numpy. I also need to verify results with propagation done for each minor body separately to check if I converted everything properly.

Initial results of @JoshPaterson propagation are awesome - entire MPCORB.dat file (~ 1M orbits) propagates for one given epoch in around 30 seconds).

[brandon-rhodes]

Wow, thank you — I look forward to reading over this tomorrow! In case you're ever interested in plotting the orbits where they look flat rather than slanted, I think the most recent Skyfield now supports instead of:

position = comets.at(t).position.au

— a maneuver like:

from framelib import ecliptic_frame
position = comets.at(t).frame_xyz(ecliptic_frame)

I only mention it because, when doing my own plots, it's always my first step so the Solar System doesn't look so lopsided. 🙂

[xmichaelx]

Thanks, it looks better:

from skyfield.api import load
from skyfield.data import mpc
from skyfield.data.mpc import _mpcorb_orbits
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from numpy import arange, linspace, rollaxis
import matplotlib.pyplot as plt
from skyfield.framelib import ecliptic_frame

ts = load.timescale()

with load.open('MPCORB.DAT') as f:
    rows = mpc.load_mpcorb_dataframe(f)
    bad_orbits = rows.semimajor_axis_au.isnull()
    rows = rows[~bad_orbits]
    rows = rows.set_index('designation', drop=False)
    comets = _mpcorb_orbits(rows.iloc[0:20], ts, GM_SUN)

    t = ts.utc(linspace(2000, 2030, 100))
    position = comets.at(t).frame_xyz(ecliptic_frame).au

    plt.rcParams['legend.fontsize'] = 10

    fig = plt.figure()
    ax = fig.gca(projection='3d')

    for i, (x, y, z) in enumerate(rollaxis(position, 1)):
        ax.plot(x, y, z, label=rows.designation[i])

    ax.legend()

    plt.show()

I'll compare results with non-vectorized version - but what would need to change in order to be able to use observe method on comets like this:

astrometric = earth.at(t).observe(sun + comets)

Current exception stack:

Traceback (most recent call last):
  File ".\sso.py", line 23, in <module>
    astrometric = earth.at(t).observe(sun + comets)
  File "d:\python-skyfield\skyfield\positionlib.py", line 583, in observe
    p, v, t, light_time = body._observe_from_bcrs(self)
  File "d:\python-skyfield\skyfield\vectorlib.py", line 106, in _observe_from_bcrs
    return _correct_for_light_travel_time(observer, self)
  File "d:\python-skyfield\skyfield\vectorlib.py", line 240, in _correct_for_light_travel_time
    tposition, tvelocity, gcrs_position, message = target._at(t)
  File "d:\python-skyfield\skyfield\vectorlib.py", line 218, in _at
    p += p2
ValueError: operands could not be broadcast together with shapes (3,) (3,20) (3,)

[brandon-rhodes]

…what would need to change in order to be able to use observe method on comets like this…

Good question. To take a single center and .observe() multiple targets, the single center’s array will need an additional dimension to match the extra dimension of the multiple targets. Immediate fixes might be either manually upgrading the dimensions of the center with something like:

p = earth.at(t)
p.shape = p.shape + (1,)

Or else improving the .observe() method so that it uses skyfields.functions._reconcile() on the center and target arrays, to add the extra missing dimension right before it starts the light-time backdating loop.

But either quick fix dodges a bigger issue that needs to be tackled: which is that this idea of a vector of comets or minor planets seems to be populating the 2nd dimension — so, the x,y,z of 20 comets has dimensions 3,20 — which conflicts with Skyfield’s design of using the 2nd dimension for time. I would have expected your sample exception to have some 100 dimensions inside, reflecting the size of the time vector, but all I see are 3 and 20 values. Did that exception come from your sample script, or a different script with a scalar rather than vector Time?

Currently Skyfield has well-founded ideas about 2 dimensions:

• 0 — Semantic: x,y,z or ra,dec,distance or whatever.
• 1 — Time: length n = length of your time t.

As the idea of “a vector of comets” expands us beyond this, we need to consider other uses of vectors, like: a vector of Earth observatories. Or a vector of other objects with which we are searching for conjunctions. A big, big question arises. Do we start prescribing for users what further dimensions always mean in Skyfield? Like:

• 2 — n centers of observation, like observatories or satellites or space probes.
• 3 — n comets or minor planets or whatever.
• 4 — n other comets, minor planets, or stars for which we are looking for conjunctions with the objects in dimension 3.

Or, since Skyfield operations all seem pretty agnostic to everything but dimensions 0 and maybe sometimes 1, do we instead leave the further dimensions 2,3,… for users to define as they wish? That’s the way I am currently leaning. We would need to write up docs, of course, for how to do it — for looking for conjunctions between comets, for example, folks would need to be told how to turn a set of positions xyz,t,n into both xyz,t,1,n and xyz,t,n,1 so than a grid of n,n mutual separations could be generated.

But a more immediate question concerns the ubiquity in current Skyfield code of the second slot always being time, except that the multi-comet and multi-asteroid code here instead uses it for n comets and asteroids. Do we add an intermediate dimension here so it comes out t,1,n and skips the “time slot” in the array dimensions? Or keep it 1,n and make this a first example where users might need to shift dimensions around if a time array gets involved later?

[xmichaelx]

But either quick fix dodges a bigger issue that needs to be tackled: which is that this idea of a vector of comets or minor planets seems to be populating the 2nd dimension — so, the x,y,z of 20 comets has dimensions 3,20 — which conflicts with Skyfield’s design of using the 2nd dimension for time. I would have expected your sample exception to have some 100 dimensions inside, reflecting the size of the time vector, but all I see are 3 and 20 values. Did that exception come from your sample script, or a different script with a scalar rather than vector Time?

Yes, sorry for that. Right now I'm testing on 20 objects, 1 epoch. And this produces above stack. I figured it would be easier to debug that full-blown many epochs/many objects case.

from skyfield.api import load
from skyfield.data import mpc
from skyfield.data.mpc import _mpcorb_orbits
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from numpy import arange, linspace, rollaxis
import matplotlib.pyplot as plt
from skyfield.framelib import ecliptic_frame

ts = load.timescale()

planets = load('de421.bsp')
earth = planets['earth']
sun = planets['sun']

with load.open('MPCORB.DAT') as f:
    rows = mpc.load_mpcorb_dataframe(f)
    bad_orbits = rows.semimajor_axis_au.isnull()
    rows = rows[~bad_orbits]
    rows = rows.set_index('designation', drop=False)
    comets = _mpcorb_orbits(rows.iloc[0:20], ts, GM_SUN)

    t = ts.utc(2021, 1, 1)
    astrometric = earth.at(t).observe(sun + comets)
    position = comets.at(t).frame_xyz(ecliptic_frame).au

    plt.rcParams['legend.fontsize'] = 10

    fig = plt.figure()
    ax = fig.gca(projection='3d')

    for i, (x, y, z) in enumerate(rollaxis(position, 1)):
        ax.plot(x, y, z, label=rows.designation[i])

    ax.legend()

    plt.show()

When 100 epochs are used stack looks like this:

Traceback (most recent call last):
  File ".\sso.py", line 24, in <module>
    astrometric = earth.at(t).observe(sun + comets)
  File "d:\python-skyfield\skyfield\positionlib.py", line 583, in observe
    p, v, t, light_time = body._observe_from_bcrs(self)
  File "d:\python-skyfield\skyfield\vectorlib.py", line 107, in _observe_from_bcrs
    return _correct_for_light_travel_time(observer, self)
  File "d:\python-skyfield\skyfield\vectorlib.py", line 246, in _correct_for_light_travel_time
    tposition, tvelocity, gcrs_position, message = target._at(t)
  File "d:\python-skyfield\skyfield\vectorlib.py", line 224, in _at
    p += p2
ValueError: operands could not be broadcast together with shapes (3,100) (3,20,100) (3,100)

So it looks like position p2 has different objects in second dimension and time is added as last dimension.

[xmichaelx]

I don't yet have a strong preference either way - although I would suspect that case where we have n objects observed will be many times more frequent than case where we have n centers of observations. For example I use sgp4 ( thank you :) ) for propagation of all the fresh TLEs I have (~20k) , convert it to ITRS manually and then check if any of them is present within given telescope FOV (by proximity to mounts target altitude and azimuth). It's hacky but works, although it would be awesome to just use skyfield for that as well.

[JoshPaterson]

But either quick fix dodges a bigger issue that needs to be tackled: which is that this idea of a vector of comets or minor planets seems to be populating the 2nd dimension — so, the x,y,z of 20 comets has dimensions 3,20 — which conflicts with Skyfield’s design of using the 2nd dimension for time.

This can be changed in keplerlb.propagate without too much trouble, so the shape would be (3, len(t), 20) instead.

[xmichaelx]

Building on what both of you've said I believe that it would be a good idea to adopt following convention:

• 0 — Semantic: x,y,z or ra,dec,distance or whatever.
• 1 — Time
• 2 — n comets or minor planets or satellites or whatever.

And leave out other proposed dimensions for now following YAGNI principle. Whenever single target is specified then last dimension can be skipped so it nicely reduces to what is available in skyfield now.

At least that's my perception of it - I'd gladly help with that but I'd need some guidance and permission from both of you :) @JoshPaterson @brandon-rhodes So, what do you think? Is support for multiple observed objects worth a try?

[brandon-rhodes]

@xmichaelx — Thanks very much for your contributions here! I look forward to getting the pull request reviewed that you inspired. (Hopefully by the end of the weekend!)

[xmichaelx]

@brandon-rhodes I'll look closely to implementation provided by @JoshPaterson and probably try to follow the same convention in order to get multiple satellites propagation using sgp4. I see that there was a related opened PR #439

@JoshPaterson awesome trick with int(, 32) :)

If you need help with test writing or anything, please let me know.