sbpy – A Python module for small-body planetary astronomy



sbpy is an astropy affiliated package that will provide tools for small-body planetary astronomers. The idea is provide a collection of well-tested and documented tools that asteroid and comet observers use on a daily basis. The goal is to improve the reproducibility of results and to make it easier for (young) researchers to try out new ideas.

The development of sbpy just started, as our NASA PDART grant just recently arrived. A lot of work to be done, but a few glimpses of its usefulness can already be seen.

The full functionality will encompass the following areas:

  • observation planning tools tailored to moving objects,
  • photometry models for resolved and unresolved observations,
  • wrappers and tools for astrometry and orbit fitting,
  • spectroscopy analysis tools and models for reflected solar light and emission from gas,
  • cometary gas and dust coma simulation and analysis tools,
  • asteroid thermal models for flux estimation and size/albedo estimation,
  • image enhancement tools for comet comae and PSF subtraction tools,
  • lightcurve and shape analysis tools, and
  • access tools for various databases for orbital and physical data, as well as ephemerides services.

If you are interested in using/testing/contributing to sbpy, check out our website and github repo for more information.


Don Quixote did it again…

Don Quixote is active again and this time the activity is observable in the optical, pointing toward the presence of dust (CBET 4502, requires a CBET subscription). More details coming soon…


In the framework of the MANOS program, I have been working on an automated data analysis pipeline for calibrated asteroid photometry and astrometry from imaging data. This pipeline, which can be applied to any kind of point source observations, is now in a robust state and available through github.

For details on the pipeline, please refer to Mommert 2017 (in press at Astronomy and Computing, arxiv).

Accessing the Gaia and Pan-STARRS catalogs using Python

The recently available Gaia and Pan-STARRS data releases are certainly interesting to the majority of astronomers for various reasons. If only there was a way to load the catalog data into Python…

Gaia DR1 catalog access using astroquery.vizier

The Gaia DR1 catalog is accessible through Vizier, which in turn can be accessed using the astroquery.vizier module, providing a comfortable astropy.table output. A simple query could look like this:

from astroquery.vizier import Vizier 
import astropy.units as u 
import astropy.coordinates as coord

def gaia_query(ra_deg, dec_deg, rad_deg, maxmag=20, 
    Query Gaia DR1 @ VizieR using astroquery.vizier
    parameters: ra_deg, dec_deg, rad_deg: RA, Dec, field 
                                          radius in degrees
                maxmag: upper limit G magnitude (optional)
                maxsources: maximum number of sources
    returns: astropy.table object
    vquery = Vizier(columns=['Source', 'RA_ICRS', 'DE_ICRS', 
                                    ("<%f" % maxmag)}, 
                    row_limit = maxsources) 
    field = coord.SkyCoord(ra=ra_deg, dec=dec_deg, 
                           unit=(u.deg, u.deg), 
    return vquery.query_region(field, 
                               width=("%fd" % rad_deg), 

# Example query
print(gaia_query(12.345, 67.89, 0.1))

Other columns are available, too. Simply add their names as provided here to the ‘columns’ list.

Pan-STARRS DR1 catalog access

The Pan-STARRS DR1 catalog resides at MAST, which unfortunately does not yet have an astroquery interface. Hence, we have to use a different approach: we download the data as an xml file and read that in using astroquery, again providing an astropy.table object:

import requests 
from import parse_single_table 
def panstarrs_query(ra_deg, dec_deg, rad_deg, mindet=1, 
    Query Pan-STARRS DR1 @ MAST
    parameters: ra_deg, dec_deg, rad_deg: RA, Dec, field 
                                          radius in degrees
                mindet: minimum number of detection (optional)
                maxsources: maximum number of sources
                server: servername
    returns: astropy.table object
    r = requests.get(server, 
    params= {'RA': ra_deg, 'DEC': dec_deg, 
             'SR': rad_deg, 'max_records': maxsources, 
             'outputformat': 'VOTable', 
             'ndetections': ('>%d' % mindet)}) 
    # write query data into local file 
    outf = open('panstarrs.xml', 'w') 
    # parse local file into astropy.table object 
    data = parse_single_table('panstarrs.xml')
    return data.to_table(use_names_over_ids=True) 
# Example query
print(panstarrs_query(12.345, 67.89, 0.1))

The file download makes this query significantly slower than comparable astroquery routines. Please note that STScI currently limits the Pan-STARRS queries on their servers to field radii smaller than 0.5 degrees.

Check Asteroid Observability with Python

Observers know the problem: there is a huge list of targets you want to observe – but what is the best time to observe them?

Use this Python script to find out:

The script requires the Python module CALLHORIZONS (see here) that can be easily installed using:

pip install callhorizons

In order to use the script, you have to modify it in three places (you can find those in the code by using a text search, looking for the word ‘user’):

  1. Target name list:
    simply create a file listing the names of your targets, one per line. You can use names, numbers, designations; make sure to use underscores in designations, e.g., 2001_AB123 instead of 2001 AB123. Each line can contain additional information on your targets, separated by a blank; this information is ignored by the script. The target name list file has to reside in the same directory as the script, unless you provide its path.
  2. Observation date range and step size:
    The script will query Horizons to obtain your targets’ ephemerides. You have to provide a start date (‘YYYY-MM-DD’; may include a UT time ‘YYYY-MM-DD HH-MM-SS’), end date, and a time step (e.g., ‘1m’ for 1 minute, ‘2h’ for 2 hours, ‘3d’ for 3 days…).
  3. Observability conditions:
    You may want to observe only targets that are brighter than a certain V magnitude, or have a maximum airmass. Enter those conditions here. See below for examples.

Run the script using:


The script will output for each target for how many time steps it is observable and write the data into two files:

  • observability.dat – observability information for each time step
  • peak_observability.dat – peak brightness information for each target over the entire period

Feel free to use, share, and modify this script. If you have questions or need help modifying the script, send me an email!

Setting Observability Conditions

You can specify you observability conditions in this line:

 observable =[((eph['V'] < 17.5) & (eph['airmass'] < 2.5))]

The conditions are given in the square brackets indexing In the given example, the target is considered observably only in those time steps when V<17.5 and when at the same time the airmass is less than 2.5. You can use all the properties specified here to setup your observability conditions.

Some examples:

observable =[((eph['V'] < 20.5) & (eph['airmass'] < 2.0) &
                       (numpy.sqrt(eph['RA_rate']**2 + 
                                   eph['DEC_rate']**2) < 0.1))]

This line also checks for the target’s total rate to be less than 0.1 arcsec/sec.

observable =[((eph['V'] < 18) & 
                      (eph['lunar_presence'] == 'dark'))]

Using this line, the script will only return those ephemerides when the Moon is below the horizon.


Spring 2016: PHY599 – Python for Scientists

Some useful links:

Timeline and scripts from the lecture:

Please note that lecture notes are no longer available here;  example scripts, however, are available on this site. For lecture notes, please refer to

  • Jan 22: some basis; data types, strings, lists
    Exercise: make Python work on your computer!
  • Jan 29: more basics; more lists, dictionaries, sets,
    Exercise: Write a code that reads this data file with Flagstaff climate data, calculates the average temperature of every month, as well as the amount of rain and snow that falls on every day of the year.  Store your results in dictionaries temp_c, rain, and snow, such that temp_c['june'] provides you the average temperature in June measured in Celsius and rain['jan 29']/snow['jan 29'] provides you the amount of rain and snow that fell on that day. Note that precipitation numbers in the data file are cumulative.
  • Feb 5: no class
  • Feb 12: advanced concepts: control flow, file i/o, modules, exceptions, classes
    Exercise: Write three different functions, each of which returns the first n numbers of the Fibonacci Sequence as a list. One of the functions has to use a for-loop, one has to use a while-loop, and the last one has to use a recursive approach (a function that calls itself over and over again).
  • Feb 19: modules, standard libraries: math, os, datetime, urllib2, subprocess
    Exercise: Write a code that grabs the NAU weather website every 5 minutes and extracts the current temperature, dew point, humidity, wind speed, and direction (hint: identify the corresponding html source code lines and extract the numbers); write the data into a file (another hint: use append when opening the file so you can interrupt your running code).
  • Feb 26: numpy: basics, arrays
    Exercise: Calculate Euler’s number e from the series of factorials (see here) in two different approaches: (1) not using any numpy functionality at all (i.e., using a simple for-loop and lists), (2) using only numpy functionality (hint: you will have to write your own factorial function, do not use the scipy one). Evaluate the series to its n-th term for n=1e2, 1e3, 1e4 and determine for each approach its runtime and the residual to numpy.e. Try to make the numpy approach as fast as possible!
  • Mar 4: more numpy: masked arrays, structured arrays, file i/o
    Exercise: Use numpy functions to read in MPCORB.DAT, the list of all known asteroids in the Solar System (here is some documentation on the structure of the file and be warned: this is a big file), into a structured array. Using as little code as possible, derive the following things:
    1. how many asteroids have q<1.3 and can be considered near-Earth asteroids? (q is the perihelion distance, which is defined as q=a(1-e)) How many are trans-Neptunian objects (q>30)?
    2. how many asteroids have accurately known orbits (uncertainty parameter U=0)? How many NEOs/TNOs?
    (Hint: only read in those columns from the file that you really need using the usecols option in genfromtxt.)
  • Mar 11: even more numpy: random, statistics, linalg, lambda function
    Exercise: Prove the Central limit theorem, by showing that for large values of lambda, the Poisson distribution approaches the normal distribution. Use numpy functions to show the convergence, e.g., using histograms.
  • Mar 18: no class (spring break)
  • Mar 25: matplotlib; Example: scatter plotting image and reference catalogs (script, image catalog, reference catalog) Exercise: implement the matching between both catalogs based on RA and Dec and try to derive the magnitude zeropoint of the image.
  • Apr 1: Example continued: match two catalogs and derive the magnitude zeropoint of an image (script, image catalog, reference catalog) Exercise:  catalog matching with the for loops will be slow for large amounts of data. Design a matching routine using numpy’s array functions (hint: you can calculate the distance of one object from the first catalog to all objects in the second catalog at one time) and then try to implement a matching routine that utilizes a kd-tree query. What are the runtimes for the invididual routines? What are the runtimes if you try to match two catalogs with coordinates created from random numbers with each other if each catalog has 1e3, 1e4, 1e5, 1e6 sources?
  • Apr 8: scipy.optimize; Example: fitting a model to data (script, data)
    Exercise: Image a binary asteroid (primary + secondary body), each of which rotates, and the secondary also orbits the primary. Derive the rotation periods of both asteroids from the data (cumulative flux from both bodies as a function of time) under the assumption that the secondary is always visible. Use the same method we used in class. Then try to implement a Lomb-Scargle Periodogram using scipy.signal.lombscargle. Do you get the same result?
  • Apr 15: scipy.interpolate; Example: interpolating 1d and 2d data (script)
  • Apr 22: jupyter, threading (threading template, example code)
  • Apr 29: no class
  • May 6: no class
  • May 13: no class


A Python Module to Query JPL HORIZONS

Working on Solar System small bodies, I make frequent use of the JPL HORIZONS system, which provides ephemerides and orbital elements at given times for any Solar System body. HORIZONS provides access through email, telnet, and a website that can be queried manually.

In order to conveniently query and investigate large amounts of data from HORIZONS, I wrote a Python module that provides access to ephemerides and orbital elements by parsing the html source code of the HORIZONS website. This method provides robust and fast access – and it’s now part of the astroquery module! Check it out!