ValueError: operands could not be broadcast together with remapped shapes

[randomwalkingg]

I am new to Skyfield and have the following question.

import numpy as np
from skyfield.api import load, wgs84

lon = np.array([0, 10, 20, 30])
ts = load.timescale()
time = ts.utc(2020, 1, 1, 0, 0, np.arange(0, 10, 1))

sites = wgs84.latlon(np.zeros_like(lon), lon, np.ones_like(lon) * 1000)
g = sites.at(time)

I got the following error:
ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (3,3,10)->(3,10,3) (3,4)->(4,3).

It is OK if do g = sites.at(time[0]). Does it mean that I cannot propagate "sites" over a time interval? And is there a way to solve my challenge without looping?

[brandon-rhodes]

Thanks for the question! This is a general NumPy issue, about which I am hoping to write up some documentation on within the next month or two. We can look at the error using simple NumPy addition and leaving Skyfield out:

import numpy as np
x = np.array([0, 10, 20, 30])
y = np.arange(0, 10, 1)
print(x + y)

The result:

ValueError: operands could not be broadcast together with shapes (4,) (10,) 

Because NumPy addition wants natively to do pairwise addition (and the same is true with multiplication and division and everything else), it runs into a problem when one array is a different length than the other.

The solution is to move the two lengths (here, 4 and 10) on to two different dimensions. NumPy is always willing to broadcast an operation when either of the lengths is 1. So here’s the solution:

import numpy as np

x = np.array([0, 10, 20, 30])
y = np.arange(0, 10, 1)

x = x[..., np.newaxis]
y = y[np.newaxis, ...]

print(x + y)

By keeping the 4 in the first axis and moving the 10 to the second axis, there's no conflict, and we get out a grid of 40×10 results:

[[ 0  1  2  3  4  5  6  7  8  9]
 [10 11 12 13 14 15 16 17 18 19]
 [20 21 22 23 24 25 26 27 28 29]
 [30 31 32 33 34 35 36 37 38 39]]

My advice would be to try the same trick with the arrays you are passing in for the time and the longitude. Fingers crossed, there's a chance that Skyfield — which is careful pretty much to only pay attention to the first axis of an array, and ignore its deeper structure — will go ahead and compute n×m positions for you.

But maybe it won't and improvements will have to be made; try out the idea and see what happens!

My guess is that keeping time as the first dimension and pushing longitude to the second dimension will work best, but please try both ways and let us know what happens. Thanks!

[randomwalkingg]

Thank you very much brandon-rhodes for your comment. Your trick really make NumPy broadcasts and does what we want. However, when applied to Skyfield, it seems that the broadcasting problem is still there.

Here are what I tried:

lon = np.array([0, 10, 20, 30])
time = ts.utc(np.arange(0, 10, 1))
lon = lon[np.newaxis, ...]
time = time[..., np.newaxis]
sites = wgs84.latlon(np.zeros_like(lon), lon, np.ones_like(lon)*1000)
g = sites.at(time)

And I got this error:

File "<ipython-input-165-5c78cee0c067>", line 1, in <module>
    g = sites.at(time)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/vectorlib.py", line 92, in at
    p, v, gcrs_position, message = self._at(t)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/toposlib.py", line 42, in _at
    RT = _T(itrs.rotation_at(t))
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/framelib.py", line 84, in rotation_at
    R = mxm(rot_z(-t.gast * tau / 24.0), t.M)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/descriptorlib.py", line 12, in __get__
    value = self.method(instance)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/timelib.py", line 751, in gast
    d_psi, _ = self._nutation_angles_radians
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/descriptorlib.py", line 12, in __get__
    value = self.method(instance)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/timelib.py", line 680, in _nutation_angles_radians
    return iau2000a_radians(self)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/nutationlib.py", line 29, in iau2000a_radians
    d_psi, d_eps = iau2000a(t.tt, fundamental_argument_terms, lunisolar_terms,
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/nutationlib.py", line 237, in iau2000a
    a = fundamental_arguments(t, fundamental_argument_terms)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/nutationlib.py", line 338, in fundamental_arguments
    a = next(fa) * t
ValueError: operands could not be broadcast together with shapes (5,1) (10,1)

If I do

lon = lon[np.newaxis, ...]
time = time[..., np.newaxis]

I get this error:

File "<ipython-input-168-5c78cee0c067>", line 1, in <module>
    g = sites.at(time)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/vectorlib.py", line 92, in at
    p, v, gcrs_position, message = self._at(t)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/toposlib.py", line 42, in _at
    RT = _T(itrs.rotation_at(t))
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/framelib.py", line 84, in rotation_at
    R = mxm(rot_z(-t.gast * tau / 24.0), t.M)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/descriptorlib.py", line 12, in __get__
    value = self.method(instance)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/timelib.py", line 751, in gast
    d_psi, _ = self._nutation_angles_radians
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/descriptorlib.py", line 12, in __get__
    value = self.method(instance)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/timelib.py", line 680, in _nutation_angles_radians
    return iau2000a_radians(self)
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/nutationlib.py", line 29, in iau2000a_radians
    d_psi, d_eps = iau2000a(t.tt, fundamental_argument_terms, lunisolar_terms,
  File "/Users/username/PythonEnv/env3/lib/python3.8/site-packages/skyfield/nutationlib.py", line 249, in iau2000a
    dpsi += dot(sarg, lunisolar_longitude_coefficients[:cutoff,1]) * t
ValueError: non-broadcastable output operand with shape (10,) doesn't match the broadcast shape (1,10)

I haven't looked at the issue #526 that you mentioned yet.

[brandon-rhodes]

Thanks for those specific exceptions! I will plan to jump in soon and figure out what Skyfield is doing in each case — hopefully I will comment here again with further information within a few days, as soon as I'm able to schedule some work on the topic. It's a direction it would be fun to see Skyfield and its documentation expand in!

[randomwalkingg]

It would be a nice feature to have in Skyfield I believe. Thanks for your helping responses.

[brandon-rhodes]

(Whoops, that commit was by accident; it named the wrong issue.)

[randomwalkingg]

@brandon-rhodes, just wondering if you have you any update on this issue?

[brandon-rhodes]

No; my recent efforts have been spent on Skyfield's estimates of where the Earth is pointing over the next few years, because its primitive ∆T curve was getting steeper and steeper every year as the Earth's rotation failed to match the rate of estimates from a decade or two ago.

But that's now fixed so I should have time within the next couple of months to return to the design/broadcasting.py file and resume experiments. While there are lots of open issues at the moment — https://github.com/skyfielders/python-skyfield/issues — I do want to give this one a bit of priority, since @JoshPaterson was generous enough to offer some pull requests with ideas about the problem.