vector calculations for apparent().altaz()

[mworion]

Hi Brandon, recently I heavily used the new panda approach for hipparchus catalogue. As I run through your examples they worked fine. Is used the same filter approach like in your documentation.

For getting an AltAz picture of the selected stars I would replace

ra, dec, distance = astrometric.radec()

straight forward with

alt, az, d = astrometric.apparent().altaz()

but this approach throws some errors:

File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/skyfield/positionlib.py", line 413, in apparent
    position_au, gcrs_position)
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/skyfield/earthlib.py", line 115, in compute_limb_angle
    coszd = dots(position_au, observer_au) / (disobj * disobs)
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/skyfield/functions.py", line 15, in dots
    return (v * u).sum(axis=0)
ValueError: operands could not be broadcast together with shapes (3,93) (3,)

From my point of view that might be related to absence of vector capabilities in apparent and altaz module. If you could confirm that, I would change my code to get a list of Stars out of the datagram like:

bright_stars = list()
for index, row in df.iterrows():
    bright_stars.append(Star.from_dataframe(row)

to get a list, which is appropriate. Any thoughts about my doing ? Many thanks in advance and merry christmas !

Michel

[brandon-rhodes]

Interesting, I didn't know that apparent() wouldn't handle this correctly! I'll see if it can be fixed.

[mworion]

Hi Brandon

out of my debugger the matrix sizes won't fit in earthlib.py method compute_limb_angle() line 115, which does: coszd = dots(position_au, observer_au) / (disobj * disobs) and even more precise: dots(position_au, observer_au) in this line throws: ValueError: operands could not be broadcast together with shapes (3,288) (3,) In my case the array observer_au is array([ 1.85504551e-05, -2.03036197e-05, 3.24741140e-05]) and positions_au is array([[ 1.80125919e+14, 1.05689414e+14, 1.98732570e+14, 4.52215070e+13, 1.51538261e+14, 1.74460494e+14, 1.11961562e+14, 1.92646015e+14, 1.07359521e+14, 9.78169398e+13, 1.35574606e+14, 1.93990307e+14, 1.59969677e+14, 9.53017087e+13, 1.39049246e+14, 1.01629986e+14, 1.78256085e+14, 1.57989068e+14, 8.03244288e+13, 1.69188053e+14, 8.53038995e+13, 1.30742238e+14, 1.60819725e+14, 1.42702733e+14, 2.08849867e+12, 1.55835323e+14, 1.12070023e+14, 1.44025913e+14, 8.49082278e+13, 1.11007542e+14, 1.06154121e+14, 8.35330204e+13, 7.80340279e+13, 3.06706045e+13, 1.02898266e+14, 9.14264777e+13, 8.02698635e+13, 1.01767711e+14, 1.00173565e+14, 4.23757243e+13, 7.69918879e+13, 4.33124625e+13, 7.09344718e+13, 6.16955487e+13, 4.68720013e+13, 3.72800146e+13, 4.49631438e+13, 3.58715542e+13, 4.63495108e+13, 4.03921174e+13, 4.02323304e+13,..... Hope this helps.

Michel

Edit: if you transpose position_au, it works fine. That's just an analysis based on matrix shapes without understanding the content ! :-(

[Bernmeister]

The following code was copied from the Skyfield manual and amended:

from skyfield.api import Star, load, Topos
from skyfield.data import hipparcos

with load.open(hipparcos.URL) as f:
    df = hipparcos.load_dataframe(f)

barnards_star = Star.from_dataframe(df.loc[87937])

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

ts = load.timescale()
t = ts.now()
astrometric = earth.at(t).observe(barnards_star)
ra, dec, distance = astrometric.radec()
print("RA:", ra)
print("Dec:", dec)

print()

boston = earth + Topos('42.3583 N', '71.0603 W')
astrometric = boston.at(t).observe(barnards_star).apparent()
az, alt, distance = astrometric.altaz()
ra, dec, distance = astrometric.radec()
print("Az:", az)
print("Alt:", alt)
print("RA:", ra)
print("Dec:", dec)

I got the following results:

RA: 17h 57m 47.48s
Dec: +04deg 44' 52.6"

Az: 00deg 35' 59.8"
Alt: 84deg 07' 08.0"
RA: 17h 57m 46.15s
Dec: +04deg 44' 49.6"

Noting the first lot of RA/DEC agree with the second...and there's no exception.

@mworion How does your code differ to mine above?

@bmxp Well spotted; now fixed.

[bmxp]

Why do you use the same args for print in

print("Az:", ra)
print("Alt:", dec)
print("RA:", ra)
print("Dec:", dec)

twice?

[mworion]

@bmxp:

your are using a single star, that work nicely, but when using a set of stars like a filtered list (in the documentation as well) like:

df = df[df['magnitude'] <= 2.5]

RaDec() works fine as well, but apparent() throws the error.

Michel

[brandon-rhodes]

So!

I've solved a more crucial bug in Skyfield and now have time to return from this. I realize that this is a new situation and I'm going to need to do a bit of paper-and-pencil thinking about it. Previously, if a position provided, say, 100 x, y, and z, positions, then it also provided 100 different times — because the model Skyfield assumes is that you get 100 positions by asking where an object like Mars is on 100 different date-times.

But this is new.

Here, you have 100 positions, but they are all for the exact same time! Instead of being for one object, as Skyfield always assumed, they are for 100 different objects.

So Skyfield, if we decide to support this situation (which would be nice?) can no longer assume that the time array is as big as the position array.

Or — or — I've got it!

Instead of rewriting all the code to now handle two different cases, we can teach the star code, when it produces 100 positions for time t, to create for them a time array with 100 copies of exactly the same time!

So maybe this doesn't need a pencil after all. Let me go take a look.

[brandon-rhodes]

Okay, I think I worked out some code! Please try:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

and let me know whether the situation — or at least the error message — improves. Thanks for letting me know about this case!

[mworion]

Brandon,

first of all, many thanks about considering this as an extension of the use of skyfield. the description you gave explains from my point of view the case. I like your pragmatic approach extending the time to the same size of the number of objects to be calculated.

The repo you mention does not exist, my guess would be to find it on the actual master branch in sky fielders GitHub location as there were the commits with your changes in. I counterchecked it, if have your changes in my environment.

Unfortunately it throws the same error message. Sizes of the matrix is still the same. Sorry to say.

Michel

[brandon-rhodes]

Ack! I should always paste URLs instead of trying to take the shortcut of hand-editing. I've fixed the link in my earlier comment.

Could you provide a script that I could try running? Since my attempt at a guess-and-fix failed, I think I'll need an actual failing case to track this further. Thanks!

[mworion]

Hi Brandon, attached the snippet for reproducing the error:

brandon.py.zip

Hope this helps. Again, many thanks for your extraordinary work! Michel

[brandon-rhodes]

Oh — I see!

The problem is that in normal calculations the observer has exactly as many positions as the observed object does: from the Earth on 10 dates we measure the distance to Jupiter on 10 other dates, and so forth. But in this case there's only 1 observer location but many object locations.

I hope there aren't a dozen new places where calculations will need re-working with this new concept of one observer but many objects observed. I guess I should try fixing this one circumstance, maybe this weekend, and see how the code behaves.

[mworion]

Exactly, The use case I have to present a list of stars for doing polar alignment of a mount. The user could select the one he would like to have (which is visible from his observing location). Michel

[bradleypallen]

Another use case is plotting stars as background to a matplotlib polar plot of a satellite pass using alt/az coordinates.

[andreww5au]

Hi Brandon, any progress on this issue?

I've just hit the same problem - if I create a Star() object with position arrays, like this:

import skyfield.api as si s = si.Star(ra=si.Angle(degrees=numpy.array([120.0, 130.0])), dec=si.Angle(degrees=numpy.array([-30.0, -40.0])))

Then I can pass it to a .observe() call to get an astrometric position, but if I call .apparent() on that astrometric position, I get an exception at line 15 in skyfield/functions.py:

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

I've been working in reverse with no problems - creating a position object using .from_altaz() with numpy arrays of thousands of Alt and Az values, and converting it with the .radec() method. That works well, and is incredibly fast compared to astropy.

If it was only a few positions, I'd just loop over them, but I'm converting existing code that created an astropy.coordinates.SkyCoord object with RA and Dec containing arrays produced by numpy.meshgrid(), with half a million positions over the whole sky.

[brandon-rhodes]

@andreww5au I just reproduced your error successfully, after fixing another bug in the vector-star case. I'll see what I can find out.

[mworion]

Hi,

first many thanks working on that. I used the actual master with your changes and tried again the provided script. It still throws errors:

Traceback (most recent call last): File "/Users/mw/PycharmProjects/MountWizzard4/gists/generate/brandon.py", line 27, in <module> alt, az, dist = observer.at(t).observe(starsDF).apparent().altaz() File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.7/site-packages/skyfield/positionlib.py", line 438, in apparent position_au, gcrs_position) File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.7/site-packages/skyfield/earthlib.py", line 116, in compute_limb_angle coszd = dots(position_au, observer_au) / (disobj * disobs) File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.7/site-packages/skyfield/functions.py", line 15, in dots return (v * u).sum(axis=0) ValueError: operands could not be broadcast together with shapes (3,15) (3,)

So for some reason I personally don't understand there might be still an open point.

Used: python 3.7.4 in virtual environment, skyfield 1.13 (but from master today) pandas 0.25.1

Michel

EDIT: Just recognized when looking on your test, that the samples from the script contain more fields that in your test: Test: ra, dec dataframe: ra shape=15, dec shape=15, ra_mas_per_year shape=15, dec_mas_per_year shape=15, parallax_mas shape=15, epoch shape=15

[andreww5au]

Hi Brandon, thanks for working on this, but I still get exactly the same error, with the latest skyfield code pulled from 'master' on github. Here's a full example showing the error.

Andrew

---

Python 3.7.3 | packaged by conda-forge | (default, Mar 27 2019, 23:01:00) [GCC 7.3.0] :: Anaconda, Inc. on linux Type "help", "copyright", "credits" or "license" for more information.

import numpy import skyfield.api as si s = si.Star(ra=si.Angle(degrees=numpy.array([120.0, 130.0])),dec=si.Angle(degrees=numpy.array([-30.0, -40.0]))) ts = si.load.timescale() [#################################] 100% deltat.data [#################################] 100% deltat.preds [#################################] 100% Leap_Second.dat t = ts.utc(2019, 10, 14) MWA_TOPO = si.Topos(longitude=(116, 40, 14.93), latitude=(-26, 42, 11.95), elevation_m=377.8) PLANETS = si.load('de421.bsp') S_MWAPOS = PLANETS['earth'] + MWA_TOPO observer = S_MWAPOS.at(t) pos = observer.observe(s) pos.radec() (<Angle 2 values from 08h 00m 00.00s to 08h 40m 00.00s>, <Angle 2 values from -40deg 00' 00.0" to -30deg 00' 00.0">, <Distance [2.06264806e+14 2.06264806e+14] au>) apos = pos.apparent() Traceback (most recent call last): File "", line 1, in File "/tmp/python-skyfield/skyfield/positionlib.py", line 438, in apparent position_au, gcrs_position) File "/tmp/python-skyfield/skyfield/earthlib.py", line 116, in compute_limb_angle coszd = dots(position_au, observer_au) / (disobj disobs) File "/tmp/python-skyfield/skyfield/functions.py", line 15, in dots return (v u).sum(axis=0) ValueError: operands could not be broadcast together with shapes (3,2) (3,)

[brandon-rhodes]

@mworion @andreww5au Interesting! The case of a Topos brings up one additional bit of math that wasn't being exercised when I tested using a plain Earth object, and an adjustment needed to be made there too.

Could you each test from master and see what exception might pop up now?

[andreww5au]

@brandon-rhodes Thanks, works perfectly now :-)

Andrew

[brandon-rhodes]

Great! I'll see about getting a release out this week, if things don't get too busy.

[mworion]

@brandon-rhodes,

Sorry for the late answer, I'm on a business trip, but got a quick test and with 1.14 is works nicely. Thanks for fixing !

Michel