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.

The module can be downloaded from github and is available through PyPi, which means you can install it using:

pip install callhorizons

Obtaining the coordinates for NEO 433 Eros as seen from Mauna Kea (observatory code 568) on a given Julian date is as easy as:

import callhorizons
eros = callhorizons.query('Eros')
eros.set_discreteepochs([2457446.177083, 2457446.182343])
print eros['RA'], eros['DEC'] 

Ephemerides and orbital elements for up to 15 discrete times, or over an arbitrary time range, can be queried for one object at a time. Data are stored in a numpy array, allowing for convenient access. Have a look at the documentation to learn more!

The module is publicly available without any warranty; feel free to use and modify it. If you have suggestions or comments, please let me know!

Observing Tools

A while back, I wrote some tools that I use for planning my observations. Now, I finally found the time to update them and turn them into publicly accessible websites. The tools are:

Observation Planner

This tool plots visibility curves indicating target elevation and airmass over the course of a whole night, which is either the current or coming night, or any night the user selects.  A number of moving and fixed targets can be displayed in the plot at a time, as shown in the example below. The purpose of this tool is to make it easier to decide which object to observe next. The black vertical represents the current time; dusk and dawn are also shown.


Pointing and Finder Tool

The idea behind this tool is to minimize the time necessary to obtain coordinates for moving targets and to create a finder chart to identify the target. Simply enter the target name and a time – or just select the current time – and the tool will provide pointing information including RA, Dec, the target’s rates, visual brightness, and airmass, as well as a finder chart indicating the target’s position at different times (see below). The tool can also be used to generate finder charts for fixed targets.



Both tools are currently set up for only a small number of observatories. If you would like to use these tools for your own observations, let me know and I can add your favorite observatory. Also, don’t hesitate to contact me if you have any other suggestions or comments!

How many Dead Comets are there?

One long-standing question is: how many asteroids – especially near-Earth asteroids (NEAs) – were comets in the past? In the classical conception, comets have surface ices and form tails and comae in the vicinity of the Sun; asteroids are simply rocks in space. In the recent years, the two classifications became less clear than that.

Comets can end up on orbits that bring them close to the Earth. We can observe them. But what happens to them after spending a while there? For some comets, like ISON, being close to the Sun is fatal; they simply disintegrate and disappear. For other comets in near-Earth space, the strong degree of insolation forces them to lose their surface ices quickly and eventually their activity stops – they are dead or dormant comets. Since comet surfaces have very low albedos, these dead comets simply look like low-albedo asteroids and are really hard to identify.

However, not all dead comets are really dead. Don Quixote is a good example, since we discovered activity in this objects that has not been detected in the last three decades. For this reason, dead comets should really be called dormant comets, as it is not clear if they ceased activity for good, or if they only sleep. The real question is: do these objects still harbor ices in their interiors? If so, this can have significant implications for the transport of water and other volatiles to Earth in the past and for the potential resource utilization in space in the future.

In order to identify dead comets that hide in the NEA population, we performed a statistical analysis in two steps. In the first step, we compared the orbits of currently active comets with those of asteroids and identified asteroids on comet-like orbits. In the second step, we identified those asteroids on comet-like orbits with comet-like (i.e., low) albedos. Those are most likely candidates for dead comets. We found a total of 23 dead comet candidates in albedo data provided by ExploreNEOs, NEOWISE, and Akari. Interestingly, only about 50% of all asteroids on comet-like orbits also have comet-like albedos.

In a second analysis, we estimated how many asteroids are likely to be of cometary origin. This question sounds easier than it is. The problem is that most asteroids are discovered with optical telescopes, which are more likely to discover asteroids with high albedos than those with low – and especially comet-like – albedos. Also, dead comets move on orbits that take them far away from the Sun, which makes them even less likely to be discovered by asteroid surveys. Hence, we have probably discovered proportionally fewer dead comets than other asteroids.

In order to resolve this discrepancy and account for this so-called ‘discovery bias’, we used a well-characterized survey that operates in the thermal-infrared, which is less prone to albedo-based discovery bias. NEOWISE is a program that uses the WISE all-sky survey data to search of all kinds of asteroids. Assuming a realistic but synthetic NEA population, we checked, how many of these synthetic NEAs would have been discovered by NEOWISE. With that information and those NEAs actually observed by NEOWISE, we can estimate how many more dead comets there are. We checked two different samples: dead comets larger than 1 km in diameter and dead comets with absolute magnitudes brighter than or equal to 21.  We find that 0.3-3.3% of the NEAs with H <= 21 and 9 (+2/-5)% of the size-limited NEA population are dead comets. These numbers are slightly lower than previously assumed, implying that fewer NEAs than previously assumed have a cometary origin.

All the details of this analysis are available in a paper that is currently in press at AJ, but already available on arxiv.