Observe from VectorSum of positions

[fgasdia]

I am trying to calculate the solar zenith angle for an array of positions. To do this I am defining a single Time and creating a VectorSum of positions by generating a Topos object with an array of lats and lons. Then I would observe().apparent().altaz(). However, I am getting a numpy error when trying to call pos.at(t).

Example:

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

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

lats = np.array([20, 30, 40])
lons = np.array([-110, -110, -110])
pos = earth + api.Topos(latitude_degrees=lats, longitude_degrees=lons)
astro = pos.at(t).observe(sun)
alt, az, dist = astro.apparent().altaz()

The problem is pos.at(t) results in TypeError: invalid data type for einsum

Is this something that could be fixed or easily supported or is the best alternative to loop over the positions?

[brandon-rhodes]

Yes, the best alternative is to loop over the positions.

I have looked into the code a bit, and even if the above error is fixed, every subsequent phase of the computation then seems to run into trouble because of the extra unexpected array dimensions encountered. Imagine, for generality, what would happen if not only your topos had an extra dimension, but that the other vectors you were composing did as well — downstream computation then can't guess what the extra dimensions represent, from what I can tell.

So I'm going to close this for the moment as infeasible, though I'll at least fix the first problem you ran into.

One big issue that would arise, API-wise: let's imagine a t with 9 dates and a Topos with 9 lat/lon pairs. The question would arise: are you asking for a single position of each Topos, on the corresponding date? Or asking for 81 positions, one for each date × position? And if the latter, which dimension would come first in the output array? I'm not sure there's a sensible way to constrain the complexity that would arise, nor to write the Python code in Skyfield so that it would handle all the possible cases correctly.

[dcajacob]

I get the same error for a different set of circumstances. In this case, I am not trying to do a vectorized operation as OP is. I find that the error is dependent on how you construct the Topos object.

target = Topos(latitude_degrees=np.degrees(lat), longitude_degrees=np.degrees(lon)) # Gives an einsum error

target = Topos('40.8939 N', '83.8917 W') # But this works

[brandon-rhodes]

@dcajacob — I haven't seen the error in that context before! If you have a moment, please open a new issue so I can tackle it. Thanks!