Support broadcast of multiple WGS84 observers against a single time and target

[andibarg]

Hi,

First of all, thanks for the great work on this repo!

I would like to observe multiple geographic positions in WGS84 from the moon. Similar to the example for observing multiple stars from earth I expect the observe function to take all positions at once (e.g. as a list). However, I get an error with the following MWE:

from skyfield.api import load, wgs84

ts = load.timescale()
t = ts.now()

planets = load('de440.bsp')
earth = planets['earth']
moon = planets['moon']

geopos = wgs84.latlon(latitude_degrees = [0,0],
                      longitude_degrees = [0,10])

moon.at(t).observe(geopos+earth)

While wgs84.latlon allows lists as inputs, the last line does not and raises the following error in "vectorlib.py" (line 219):

ValueError: operands could not be broadcast together with shapes (3,) (3,2)

Is this a bug? I cannot find an easy fix but it seems that the shapes should be compatible somehow.

Cheers

[brandon-rhodes]

Good question! It is not a bug, strictly speaking, because Skyfield only promises to be able to "broadcast" (to use the NumPy term for a single operation that's made to span a whole array) along the time axis, providing support for multiple times. But it does not provide any general mechanisms for broadcasting across multiple planets, or observatories, or Earth satellites. On a few occasions I have wished that it could easily be expanded to do so, but I have never worked out a way to convince the underlying NumPy library to play along.

In this case, the problem happens even before Skyfield can really get started doing the light-time back-dating of the Earth positions. Even just trying to compute a multi-observatory position fails! So replacing the last line of your script with this:

(earth + geopos).at(t)

— is enough to trigger the exception.

And I don't see any way to — to —

twenty minutes of fiddling ensue

— well, drat.

This almost ended so simply. I almost said, "NumPy doesn't want to cooperate, so what can we do," and so work with you to at least make an acceptably efficient loop for you that would maybe not re-compute the same Earth position n times for n observatories, but would otherwise solve the problem with an old-fashioned for loop.

But then I vaguely remembered having written a routine that adjusts NumPy arrays so that they broadcast correctly, even when one of them doesn't have enough dimensions. Which is what happens in this case, where there's only one Earth, but 2 observatories.

So I went looking and found the routine. It's named _reconcile(). So I started sticking it in everywhere your code raised an exception, and, to my surprise, it didn't need to be added in one hundred different places, but only in three places, and all in the same file. This might not render well, but here's what I tried:

diff --git skyfield/vectorlib.py skyfield/vectorlib.py
index a571eb5..2986f8e 100644
--- skyfield/vectorlib.py
+++ skyfield/vectorlib.py
@@ -5,7 +5,7 @@ from numpy import max
 from .constants import C_AUDAY
 from .descriptorlib import reify
 from .errors import DeprecationError
-from .functions import length_of
+from .functions import length_of, _reconcile
 from .positionlib import build_position
 from .timelib import Time
 
@@ -215,8 +215,10 @@ class VectorSum(VectorFunction):
             p2, v2, _, message = vf._at(t)
             if vf.center == 399:
                 gcrs_position = -p
-            p += p2
-            v += v2
+            p, p2 = _reconcile(p, p2)
+            p = p + p2
+            v, v2 = _reconcile(v, v2)
+            v = v + v2
         if vfs[0].center == 0 and vf.center == 399:
             gcrs_position = p2
         return p, v, gcrs_position, message
@@ -241,6 +243,7 @@ def _correct_for_light_travel_time(observer, target):
 
     tposition, tvelocity, gcrs_position, message = target._at(t)
 
+    tposition, cposition = _reconcile(tposition, cposition)
     distance = length_of(tposition - cposition)
     light_time0 = 0.0
     for i in range(10):
@@ -259,6 +262,7 @@ def _correct_for_light_travel_time(observer, target):
         light_time0 = light_time
     else:
         raise ValueError('light-travel time failed to converge')
+    tvelocity, cvelocity = _reconcile(tvelocity, cvelocity)
     return tposition - cposition, tvelocity - cvelocity, t, light_time
 
 def _jpl_name(target):

And here's your script, expanded to also produce the observatory positions individually, to double-check them:

from skyfield.api import load, wgs84

ts = load.timescale()
t = ts.now()

planets = load('de440.bsp')
earth = planets['earth']
moon = planets['moon']

geopos = wgs84.latlon(latitude_degrees = [0,0],
                      longitude_degrees = [0,10])

pos = moon.at(t).observe(geopos+earth)
print(pos.xyz)

pos1 = moon.at(t).observe(earth + wgs84.latlon(
    latitude_degrees = 0,
    longitude_degrees = 0,
))
print(pos1.xyz)

pos2 = moon.at(t).observe(earth + wgs84.latlon(
    latitude_degrees = 0,
    longitude_degrees = 10,
))
print(pos2.xyz)

And the result:

[[-0.00241964 -0.00242047]
 [-0.00053748 -0.00054486]
 [-0.00030687 -0.00030687]] au
[-0.00241964 -0.00053748 -0.00030687] au
[-0.00242047 -0.00054486 -0.00030687] au

So it looks like your case, which I'd never thought of before, might be able to be supported with just a few lines of new code, that would wind up as nearly no-op's under normal circumstances since the vectors wouldn't need adjustment, but that in your case would do the bit of extra work to make NumPy able to broadcast the n observatories.

I'll have to think a bit to make sure I'm not digging a hole here that Skyfield later won't be able to climb out of, but offhand I don't see any obvious reason not to add those lines to the next release. I'll give myself a couple of days to mull over it.

[andibarg]

Thank you so much for your quick response. I am happy to see that there could be support with just a few extra lines of code. Of course, it is important that skyfield does not break elsewhere because of this.

By chance I noticed that now adding .apparent() raises a similar error as before, but this time it is because of the numpy.einsum function in relativity.py (line 113):

ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (3,2)->(2,3) (3,3)->(3,3)

I tried to make it work but could not, as I am not too familiar with Einstein summation. Do you happen to know?

[andibarg]

I found that the error is due to observer_position in relativity.py (line 113) having the wrong shape. Instead of (3,) its shape is (3,3) if we use multiple lats/longs. However, it looks like the diagonal of the 3-by-3 array is the one we actually want. So adding the following lines in the beginning of the light_time_difference function in relativity.py solves the error and gives the correct results for the example:

if observer_position.shape == (3,3):
    observer_position = np.diagonal(observer_position)

Unfortunately, there is still an issue if I want to observe the geographic positions with an EarthSatellite. In order to make it work in this case, I need to also modify the file earthlib.py in line 94:

-    coszd = dots(position_au, observer_au) / (disobj * disobs)
+    coszd = np.dot(observer_au, position_au) / (disobj * disobs)

I could not make it work with the dots function and needed to use the numpy one.

[andibarg]

I implemented the changes discussed here in pull request #1041. Hope this is ok. I would love to see this functionality as part of skyfield.

[brandon-rhodes]

I implemented the changes discussed here in pull request #1041. Hope this is ok. I would love to see this functionality as part of skyfield.

No worries — I don't mind having my example code turned into a PR!

While most of my focus this month has been on trying to finally land an obscure fix to the sunrise routine so that I can finally put that effort to rest before starting to make other changes to the code, I have been weighing this issue in the back of my mind. The only hesitation I have is that it's a mismatch for how NumPy works. In NumPy, the user manually runs .reshape() to make their arrays broadcast correctly, and it might be nice for Skyfield to match that explicit behavior, rather than trying to magically guess what the caller wants.

And it's easy enough to write a .reshape() method for the position class.

The problem comes when your code wants to .observe(geopos+earth) because observe() starts asking the geopos+earth over and over again for positions that are earlier and earlier in time until it figures out where the position was when the light left it that's now arriving at the observer. And earth returns a single position for time t, but geopos returns two, and so there's a mismatch right down inside the + operation where the caller doesn't have time to do a reshape(), because the caller never gets their hands on those intermediate results.

So I need to do some investigation of what it would look like if there were a way for the user to ask for the earth positions to get reshaped. If there's an easy way, then I will probably opt for the more explicit design where users do reshaping themselves. If not — then maybe we'll try the magic lines I wrote.

One problem with trying to do it by magic is: what if there are several axes? Your case is simple:

(3,) earth
(3, 2) observatories

So simply extending the earth shape to (3, 1) makes the arrays able to broadcast.

But I recall discussions where folks have wanted to do things like n satellites from m observatories. In that case we have:

(3, n) satellites
(3, m) observatories

And what needs to happen is reshaping that moves those axes away from each other:

(3, 1, n) satellites
(3, m, 1) observatories

And I'm not sure that's something that Skyfield should be guessing about? But probably under the control of the user? Because, if it happens that n == m, we can't tell if:

• They want an n × n array of every satellite from every observatory.
• Or, each observatory is being paired with exactly one corresponding satellite, and they want no reshaping done.

So if we can work out a reshaping approach, then it will be the more flexible going forward. And will make it much easier to determine what's going on in people's code, if it doesn't invisibly create shapes that will then come out in error messages without the user knowing where they're coming from.

Anyway, that's the state of my thoughts at the moment!