Automatically create finding charts for known asteroids and comets with python

rock_finding_chart.md

This code will display finding charts for moving bodies (e.g., comets, asteroids) using either DS9 or Matplotlib. The code will:

• query JPL Horizons to determine the R.A. and Dec. of the target at the given time; default = Time.now()
• download a POSS2 Red image from the STScI DSS server
• display the image in either DS9 or Matplotlib
• mark the expected location of the object with a green circle
• mark the past/future locations of the object with purple/red circles in 15 minute increments
• print the R.A. and Dec. coordinates at the requested time
• check and print a warning if the Moon is within 30 degrees of the target

All you need to do is specify the name of the target in a format that Horizons can understand (e.g., target_id = 'C/2019 Q4').

rock_finding_chart.py

import numpy as np
from astropy.time import Time
from astropy.coordinates import SkyCoord, get_moon
from urllib.parse import urlencode
from astropy.io import fits
from astropy.stats import sigma_clip
from astropy import wcs
import astropy.units as u
from astroquery.jplhorizons import Horizons
import matplotlib.pyplot as plt

def get_ephemerides(target_id, id_type, time=None, location=None):
    """
    wrapper to query Horizons and return the object ephemerides.
    
    With `id_type = 'smallbody'`, you can specify the object by name
    e.g.,  `target_id = 'C/2019 Q4'` or `target_id = '2000 WK63'`
    
    With `id_type = 'id'`, you can specify the object by its ID number
    e.g., `target_id = '90000392'`
    
    `time` -- (optional) should be specified as an astropy.time.Time object
    
    `location` -- (optional) should be an astropy.coordinates.EarthLocation
    """
    if time is None:
        # use current time by default
        time = Time.now()   
    obj = Horizons(id=target_id, id_type=id_type, location=location, epochs=time.jd)
    return obj.ephemerides()

def get_dss(pos, width=15.0, height=15.0, survey='poss2ukstu_red'):
    """
    retrieve a FITS image from the STScI Digitized Sky Survey server
    
    `pos` is a SkyCoord specifying the center of the output image
    
    `width` and `height` give the dimensions of the output image in arcminutes
    
    Use `survey` to specify the survey name and band
    """
    # build the URL to download the image
    url = "http://archive.stsci.edu/cgi-bin/dss_search?"
    options = {}
    options['v'] = survey
    ra, dec = pos.to_string('hmsdms', sep=':').split()
    options['r'] = ra
    options['d'] = dec
    options['e'] = 'J2000'
    options['h'] = height
    options['w'] = width
    options['f'] = 'fits'
    query = urlencode(options)
    return fits.open(url + query)

def make_chart(target_id, id_type, time=None, useDS9=True):
    """
    create a simple finding chart for a moving object.
    
    `target_id` and `id_type` are used to retrieve the objects epherimis from Horizons
    
    The position of the target at the specified `time` (defaults to astropy.time.Time.now())
    is marked on an image of the field with a green circle. Circles marking the motion of the
    target in 15 minute intervals are added. Red circles mark future positions and purple
    circles mark past positions. The target name, observation epoch (for the green circle), 
    and coordinates at the specified `time` are printed. 
    
    If `useDS9` is False, the chart will be drawn using matplotlib instead of DS9.
    """
    
    # default to the current time
    if time is None:
        time = Time.now()

    # get the position of the target at this time
    eph = get_ephemerides(target_id, id_type, time)
    target_name = eph['targetname'][0]
    pos = SkyCoord(eph['RA'], eph['DEC'], unit='deg')[0]

    # tell the user how fast the target is moving
    message = "{} is moving at: {:.2f} \"/minute in RA; {:.2f} \"/minute in Dec;"
    print(message.format(target_name, eph['RA_rate'][0] / 60, eph['DEC_rate'][0] / 60))

    # download the image 
    poss2 = get_dss(pos) #, survey='poss1_red')
    imwcs = wcs.WCS(poss2[0].header)

    if useDS9:
        # show in DS9
        import pyds9
        ds9 = pyds9.DS9()
        ds9.set_fits(poss2)
        ds9.set('scale zscale')
        ds9.set('zoom to fit')

        # Label the chart
        ny, nx = poss2[0].data.shape
        ds9.set('regions delete all' )
        label = 'image; text {} {} # text={{{}}} font={{courier 14 bold}}'
        ds9.set('regions', label.format(nx/2, 0.95*ny, target_name))
        ds9.set('regions', label.format(nx/2, 0.90*ny, pos.to_string('hmsdms')))
        ds9.set('regions', label.format(nx/2, 0.85*ny, '(at {} UT)'.format(time.iso)))

        # mark the expected location of the target at this time
        x, y = imwcs.wcs_world2pix(pos.ra.deg, pos.dec.deg, 1)
        ds9.set('regions', 'image; circle {} {} {} # color=green'.format(x, y, 10))
    else:
        # display the image using Matplotlib
        image = poss2[0].data
        bg = sigma_clip(image)
        vmin = bg.mean() - bg.std()
        vmax = bg.mean() + 5 * bg.std()
        plt.figure(figsize=(15,15))
        plt.subplot(projection=imwcs)
        plt.imshow(image, origin='lower', cmap='gray', vmin=vmin, vmax=vmax)
        plt.xlabel('R.A.')
        plt.ylabel('Dec.')

        # Label the chart
        opts = {}
        opts['horizontalalignment'] = 'center'
        opts['verticalalignment'] = 'center'
        opts['transform'] = plt.gca().transAxes
        opts['color'] = 'green'
        opts['fontsize'] = 18
        opts['fontweight'] = 'bold'
        opts['bbox'] = dict(facecolor='black', alpha=0.75)        
        plt.text(0.5, 0.95, target_name, **opts)
        plt.text(0.5, 0.90, pos.to_string('hmsdms'), **opts)
        plt.text(0.5, 0.85, '(at {} UT)'.format(time.iso), **opts)
        
        # mark the expected location of the target at this time
        x, y = imwcs.wcs_world2pix(pos.ra.deg, pos.dec.deg, 0)
        plt.plot(x, y, 'go', mfc='None', ms=20)

    # Moon proximity warning
    moon = get_moon(time)
    sep = moon.separation(pos)
    if sep < (30 * u.deg):
        warn = "Moon is {:.1f} away!".format(sep)
        print(warn)
          
        if useDS9:
            label = 'image; text {} {} # text={{{}}} font={{courier 14 bold}} color={{red}}'
            ds9.set('regions', label.format(nx/2, 0.80*ny, warn))
        else:
            opts['color'] = 'red'
            plt.text(0.5, 0.80, warn.format(time.iso), **opts)

    # mark the previous/future locations of the target
    tdeltas = np.arange(-60, 61, 15)
    w = tdeltas != 0
    times = time + tdeltas[w] * u.min
    ephs = get_ephemerides(target_id, id_type, times)
    for i, thiseph in enumerate(ephs):
        if tdeltas[i] < 0:
            color = 'purple'
        else:
            color = 'red' 
        pos = SkyCoord(thiseph['RA'], thiseph['DEC'], unit='deg')
          
        if useDS9:
            x, y = imwcs.wcs_world2pix(pos.ra.deg, pos.dec.deg, 1)
            ds9.set('regions', 'image; circle {} {} {} # color={}'.format(x, y, 10, color))    
        else:
            x, y = imwcs.wcs_world2pix(pos.ra.deg, pos.dec.deg, 1)
            plt.plot(x, y, marker='o', color=color, mfc='None', ms=20)


## Examples

# set the target name for Horizons
# e.g.:
#   target_id = 'C/2019 Q4'
#   target_id = '2000 WK63'
#   target_id = '4450'
target_id = '2020 AV2'
id_type = 'smallbody'

# *OR* set the target id
# e.g.:
#   target_id = '90000392'
#   id_type='id'

# use DS9 to show the chart
make_chart(target_id, id_type)

# use matplotlib to display the chart
make_chart(target_id, id_type, useDS9=False)