skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

Exception when computing comet orbit #488

Closed Bernmeister closed 3 years ago

Bernmeister commented 3 years ago

When running the code below, I get the following exception on a bunch of comets:

The truth value of a Series is ambiguous. Use a.empty, a.bool(), a.item(), a.any() or a.all().

The code can either load the comet data from the Minor Planet Center OR use a locally saved copy (to avoid hammering the MPC server).

The test code:

#!/usr/bin/env python3

# Download and save https://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft00Cmt.txt

import datetime, io, skyfield.api, skyfield.constants, skyfield.data.mpc

utcNow = datetime.datetime.strptime( "2020-11-24", "%Y-%m-%d" )
latitude = -33
longitude = 151
elevation = 100

print( "Skyfield:", skyfield.__version__ )

timeScale = skyfield.api.load.timescale( builtin = True )
t = timeScale.utc( utcNow.year, utcNow.month, utcNow.day, utcNow.hour, utcNow.minute, utcNow.second )
ephemeris = skyfield.api.load( "de421.bsp" )
sun = ephemeris[ "sun" ]
earth = ephemeris[ "earth" ]
topos = skyfield.api.Topos( latitude_degrees = latitude, longitude_degrees = longitude, elevation_m = elevation )
alt, az, earthSunDistance = ( earth + topos ).at( t ).observe( sun ).apparent().altaz()

download = False
if download:
    with skyfield.api.load.open( skyfield.data.mpc.COMET_URL ) as f:
        dataframe = skyfield.data.mpc.load_comets_dataframe( f ).set_index( "designation", drop = False )

else:
    with open( "Soft00Cmt.txt" ) as f:
        orbitalElementData = f.readlines()

    with io.BytesIO() as f:
        for orbitalElement in orbitalElementData:
            f.write( ( orbitalElement + '\n' ).encode() )

        f.seek( 0 )
        dataframe = skyfield.data.mpc.load_comets_dataframe( f ).set_index( "designation", drop = False )

for name, row in dataframe.iterrows():
    try:
        body = sun + skyfield.data.mpc.comet_orbit( dataframe.loc[ name ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )
        ra, dec, earthBodyDistance = ( earth + topos ).at( t ).observe( body ).radec()
        ra, dec, sunBodyDistance = sun.at( t ).observe( body ).radec()
    except Exception as e:
        print( name )
        print( e )

I have also attached the data file (in case there is something odd in today's data and it disappears tomorrow).

Soft00Cmt.txt

brandon-rhodes commented 3 years ago

@Bernmeister — You have run into a subtle behavior of Pandas. If .loc[name] asks about a name for which only one row exists, you get back that single row as a Series. But if the name exists several times in the file, you get back a DataFrame with all n of the matching rows! Yes, a grand source of bugs: depending on the data, it can return either of two different types, so that the code that follows must be prepared for two different data types.

I confirmed the problem by adding print statements until the behavior was clear. Here's what I did to your loop, in case you need to do further investigations like this in the future:

for name, row in dataframe.iterrows():
    print(row)
    print(row.eccentricity)
    print(dataframe.loc[ name ])
    print(dataframe.loc[ name ].eccentricity)
    if dataframe.loc[ name ].eccentricity == 1.0:
        print('e == 1.0')

Catching the exception, by the way, was hiding the traceback, which is necessary to seeing where the code is failing; so as my first step I removed the try…except maneuver to learn more about what was going on.

I would suggest simply using the row value you already have instead of indexing back into the dataframe with .loc[], which goes to the expense of building the row object all over again.

Try .comet_orbit(row) and see if your script's behavior improves!

Bernmeister commented 3 years ago

Well that worked! Still trying to get my head around numpy (and frankly only heard about it once I started using Skyfield). Thanks again!