Asteroid Detection Module

Lecturer: Ashish Mahabal
Jupyter Notebook Authors: Quanzhi Ye, Ashish Mahabel, Dmitry Duev, & Cameron Hummels

This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html

Objective

Compute the detection rate of fast-moving asteroids, asteroids that pass within ~20 lunar distances from Earth, and leave trails/streaks on astronomical images.

Key steps

  • Figure out which known asteroids will show up as streaks in ZTF images, and
  • Examine the streak catalog to see how many of them are detected.

Required dependencies

See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with pip install <module>. The external astromatic packages are easiest installed using package managers (e.g., rpm, apt-get).

Python modules

  • python 3
  • astropy
  • numpy
  • astroquery
  • matplotlib

External packages

None

In [1]:
import requests, numpy as np, glob, matplotlib.pyplot as plt, os, tarfile
from astroquery.jplhorizons import Horizons
from astropy.time import Time
from IPython.display import Image, display
import pdb
import warnings

Step 1: fast-moving asteroids that fall on ZTF images

We will loop over the catalog of known NEAs (Near-Earth Asteroids) and make a series of queries to JPL's Horizon service and IPAC's IRSA service, to determine if and when they will be visible at Palomar (where ZTF runs from).

Let's download the catalog of known NEAs from the Minor Planet Center (MPC).

In [2]:
os.chdir('data') # All relevant data are in data directory
print('retrieving MPCORB...')
# mpcorb = 'https://www.minorplanetcenter.net/iau/MPCORB/NEA.txt'
# r = requests.get(mpcorb)
# mpcorb = r.text.split("\n")
mpcorb = open('NEA.txt','r').read().split("\n")
print('MPCORB file retrieved.')
retrieving MPCORB...
MPCORB file retrieved.

Now let us loop over the MPC catalog to check the visibility for every NEA. We will do it for NEA 2018 CL (the first NEA found by ZTF), and it's up to you to write a loop for it.

The format of the MPC catalog is described here, it's best for machines but is somewhat human readable. First, we create an instance for the Horizons query. We only need to tell Horizons the name of the object we want to query, and the time window we want to query, and let Horizons do the heavy-lifting. Oh, and the three-letter-code for ZTF is I41, as described here. If the query returns zero entries, it means that the object is not observable at Palomar at the suggested time (interval).

First, let us define a time window we would like to query for. Our default is 20180205, the day ZTF science observations started:

In [3]:
start_date = '20180205'
end_date = '20180206'

Then let's do all the tasks described above. For now, we do our test with 2018 CL (the first NEA found by ZTF). (in your local copy of astropy - astroquery/jplhorizons/init.py - you may have to make a change. Replace http below with https: 'http://ssd.jpl.nasa.gov/horizons_batch.cgi',)

In [4]:
def get_mpc_primary(inp):
    """ extract primary ID from an MPC-1 string parameter
    ---------
    input: MPC-1 data, strip of column 166-190
    """

    if inp[0:1] == '(':
        # (433) Eros
        return inp[inp.find("(")+1:inp.find(")")]
    else:
        # 2018 CL
        return inp

for i, obj in enumerate(mpcorb):
    obj_name = get_mpc_primary(obj[166:190].strip())
    if not obj_name == '2018 CL':
        continue
    else:
        break

obj = Horizons(id=obj_name, location='I41', \
                    epochs={'start': '{}-{}-{}'.format(start_date[0:4], start_date[4:6], start_date[6:8]), \
                   'stop': '{}-{}-{}'.format(end_date[0:4], end_date[4:6], end_date[6:8]), \
                   'step': '1h'})

eph = obj.ephemerides(skip_daylight=True, airmass_lessthan=5.0)
if len(eph) <= 0:
    print('object {} was not observable at Palomar'.format(obj_name))

We also want to calculate the motion of the asteroid to see if it will show up as a streak in our images.

In [5]:
rate = np.sqrt(eph['RA_rate']**2+eph['DEC_rate']**2)

Now, since we just discovered the power of JPL Horizons... why don't we take a short detour and do something fun? Let's say, we want to see the predicted trajectory of 2018 CL...

In [6]:
fig = plt.figure()
ax = fig.add_subplot(111)

ax.plot(eph['RA'], eph['DEC'], 'o-')
ax.set_title('trajectory of 2018 CL on from %s to %s' % (start_date, end_date))
ax.set_xlabel('RA (deg)')
ax.set_ylabel('DEC (deg)')

for i, xy in enumerate(zip(eph['RA'], eph['DEC'])):
    ax.annotate('%s' % eph['datetime_str'][i], xy=xy, textcoords='data')

plt.show()

What about the change of azimuth and elevation angle?

In [7]:
fig = plt.figure()
ax = fig.add_subplot(111)

ax.plot(eph['AZ'], eph['EL'], 'o-')
ax.set_title('azimuth/elevation angle of 2018 CL on from %s to %s' % (start_date, end_date))
ax.set_xlabel('Azimuth (deg)')
ax.set_ylabel('Elevation (deg)')

for i, xy in enumerate(zip(eph['AZ'], eph['EL'])):
    ax.annotate('%s' % eph['datetime_str'][i], xy=xy, textcoords='data')

plt.show()

The complete lists of the values supported by astroquery.jplhorizons can be found here. Play with it and be creative! For example, can you plot the change of RA/DEC rate over time? And what about the change of predicted V magnitude?

Now let's come back to our main journey. How fast does the asteroid need to be in order to show up as a streak? Let's take anything longer than 5 full-width-half-maximum (FWHM) as a streak. We know the exposure time of ZTF is 30 seconds, and the typical FWHM is 2''. Therefore, the minimum rate of motion is not too difficult to calculate:

In [8]:
min_rate = 2*5/30

print('minimum rate of motion for a streak: %.2f arcsec/sec.' %  min_rate)
minimum rate of motion for a streak: 0.33 arcsec/sec.

If the object will be a fast-mover at some point within the time window being queried, we will use IPAC's Moving Object Search Tool (MOST) to get a list of the images that the object is a fast-mover on them.

In [9]:
if max(rate) > min_rate:
    params = {'catalog': 'ztf', 'input_type': 'name_input', 'obj_name': '{}'.format(obj_name), \
              'obs_begin': '{}-{}-{}'.format(start_date[0:4], start_date[4:6], start_date[6:8]), \
              'obs_end': '{}-{}-{}'.format(end_date[0:4], end_date[4:6], end_date[6:8]), \
              'output_mode': 'Brief'}
    irsa_return = requests.get('https://irsa.ipac.caltech.edu/cgi-bin/MOST/nph-most', params=params)

...and MOST did give us something, right?

In [10]:
print(irsa_return.text)
\output_url = "https://irsa.ipac.caltech.edu:443/workspace/TMP_17w7Dr_2536/MOST/pid2536"
\catalog = "ztf"
\user_input_begin_time = "2018-02-05"
\user_input_end_time = "2018-02-06"
\object_name = "(2018 CL)"
\element_epoch = "58154.0000   (2018-Feb-05 00:00:00.0000)"
\eccentricity = " 0.242496961324341"
\inclination = "11.987721695251381"
\argument_perihelion = "142.423969683856910"
\ascending_node = "136.338359646885010"
\semimajor_axis  =  " 0.851979508669249"
\mean_anomaly = "236.573041656847607"
\magnitude_parameters = "25.50  0.00" 
\job_time_stamp = "Mon Apr  1 22:42:14 2019"
\description = ZTF processed science image
\identifier  = ivo://irsa.ipac/ztf.sci
\notes       = ZTF science images
\datatype    = catalog
\archive     = IRSA
\set         = ZTF
\ generated at Mon Apr  1 08:00:01 PDT 2019
\                                                                               
\fixlen = T
\ 
\ ra (deg) 
\ ___ Right ascension of image center
\ dec (deg) 
\ ___ Declination of image center
\ filefracday 
\ ___ Integer (YYYYMMDDdddddd) used in product file names
\ field 
\ ___ ZTF field number
\ ccdid 
\ ___ CCD number (1..16)
\ qid 
\ ___ Quadrant ID (1..4)
\ rcid 
\ ___ Readout-channel ID (0..63)
\ fid 
\ ___ Filter ID
\ filtercode 
\ ___ Filter code (abbreviated name)
\ pid 
\ ___ Science product ID
\ nid 
\ ___ Day/night ID
\ expid 
\ ___ Exposure ID
\ itid 
\ ___ Image type ID
\ imgtypecode 
\ ___ Single letter image type code
\ obsdate 
\ ___ Observation UT date/time
\ obsjd (d) 
\ ___ Observation Julian day
\ ra1 (deg) 
\ ___ Right ascension of first image corner
\ dec1 (deg) 
\ ___ Declination of first image corner
\ ra2 (deg) 
\ ___ Right ascension of second image corner
\ dec2 (deg) 
\ ___ Declination of second image corner
\ ra3 (deg) 
\ ___ Right ascension of third image corner
\ dec3 (deg) 
\ ___ Declination of third image corner
\ ra4 (deg) 
\ ___ Right ascension of fourth image corner
\ dec4 (deg) 
\ ___ Declination of fourth image corner
\ 
|    ra_obj|   dec_obj|sun_dist|geo_dist|  dist_ctr|   phase|  vmag|match|                ra|               dec|   filefracday|       field| ccdid|  qid| rcid|  fid| filtercode|                  pid|   nid|       expid| itid| imgtypecode|                       obsdate|             obsjd|               ra1|              dec1|               ra2|              dec2|               ra3|              dec3|               ra4|              dec4|
|    double|    double|  double|  double|    double|  double|double|  int|            double|            double|          long|         int|   int|  int|  int|  int|       char|                 long|   int|         int|  int|        char|                          char|            double|            double|            double|            double|            double|            double|            double|            double|            double|
 130.243140  10.591412   0.9936   0.0077     0.3695   9.6662  15.69     1   129.946762820000    10.364530940000 20180205231528          517      6     4    23     2          zr          400231532315    400     40023153     1            o            2018-02-05 05:33:25  2458154.731539400   130.387829180000    10.796160510000   129.508055050000    10.798784690000   129.506884090000     9.932328820000   130.383994010000     9.929949140000 
 129.946506  11.101031   0.9935   0.0077     0.1305   9.6479  15.67     1   129.949891860000    11.231928510000 20180205251296          517      6     1    20     2          zr          400251292015    400     40025129     1            o            2018-02-05 06:01:52  2458154.751296300   130.392207100000    11.663277090000   129.509976810000    11.666108690000   129.508681130000    10.799672070000   130.388448090000    10.797038270000 
 129.625980  11.645721   0.9934   0.0076     0.5217   9.6703  15.65     1   129.950695400000    11.232704280000 20180205272106          517      6     1    20     2          zr          400272102015    400     40027210     1            o            2018-02-05 06:31:50  2458154.772106500   130.392985830000    11.664048400000   129.510776420000    11.666846170000   129.509495230000    10.800441100000   130.389245740000    10.797862500000 
 129.618211  11.658854   0.9934   0.0076     0.3954   9.6713  15.65     1   129.648394110000    12.053413850000 20180205272593         1563      2     3     6     2          zr          400272600615    400     40027260     1            o            2018-02-05 06:32:33  2458154.772604200   130.091191430000    12.485962890000   129.206230010000    12.486622890000   129.207327970000    11.620350100000   130.089107560000    11.619537940000 
 129.304417  12.186967   0.9934   0.0075     0.2490   9.7350  15.64     1   129.070327660000    12.285688820000 20180205292488          517     10     3    38     2          zr          400292483815    400     40029248     1            o            2018-02-05 07:01:11  2458154.792488400   129.513283440000    12.718624570000   128.627359850000    12.718548800000   128.629102050000    11.852312030000   129.511895640000    11.852222940000 
 129.296470  12.200286   0.9934   0.0075     0.3748   9.7371  15.64     1   129.649259050000    12.054213010000 20180205292986         1563      2     3     6     2          zr          400292980615    400     40029298     1            o            2018-02-05 07:01:54  2458154.792986100   130.092073220000    12.486767920000   129.207114190000    12.487391270000   129.208189100000    11.621136410000   130.089970620000    11.620309610000 
 128.976880  12.733908   0.9933   0.0075     0.4280   9.8423  15.62     1   129.071104030000    13.152379000000 20180205312789          517     10     2    37     2          zr          400312803715    400     40031280     1            o            2018-02-05 07:30:26  2458154.812800900   129.515614710000    13.585218500000   128.626757470000    13.584987180000   128.628208840000    12.718695140000   129.514139550000    12.718795740000 

Here is a loop over the query result from MOST to find out which images will contain the target as a fast-mover. Note that we want the predicted magnitude of the object to be <20 since that's the typical depth of ZTF.

In [11]:
limiting_magnitude = 20

We also make another query to Horizons to ensure the object is really a fast-mover in these dates. We do this as a double-check since Horizons produces most precise ephemerides for NEAs. Horizons also tell us what the expected positional uncertainty will be. Obviously, we don't want the positional uncertainty to be too big -- say, over 20".

In [12]:
max_unc = 20

And then we print out all the metadata that are potentially useful. We will just return the first entry (or it will be a very long entry!)

In [13]:
for ztf_entry in irsa_return.text.split("\n"):
    try:
        ztf_entry_vmag = float(ztf_entry[62:67])
    except:
        continue

    if ztf_entry_vmag < limiting_magnitude:
        ztf_entry_jd = float(ztf_entry[270:287])
        ztf_entry_pid = int(ztf_entry[185:199].strip())
        ztf_entry_fracday = int(ztf_entry[120:126].strip())
        ztf_entry_filefracday = int(ztf_entry[112:126].strip())
        ztf_entry_field = int(ztf_entry[135:139].strip())
        ztf_entry_filtercode = str(ztf_entry[173:177]).strip()
        ztf_entry_ccdid = int(ztf_entry[144:146].strip())
        ztf_entry_imgtypecode = str(ztf_entry[235:238]).strip()
        ztf_entry_qid = int(ztf_entry[150:152].strip())
        ztf_entry_sundist = float(ztf_entry[24:31])
        ztf_entry_geodist = float(ztf_entry[34:40])
        ztf_entry_dateUT = str(ztf_entry[112:120]).strip()
        
        obj_this = Horizons(id=obj_name, location='I41', epochs=ztf_entry_jd)
        eph = obj_this.ephemerides()
        ztf_entry_rate = np.sqrt(eph['RA_rate'][0]**2+eph['DEC_rate'][0]**2)

        if ztf_entry_rate > min_rate:
            ztf_entry_unc = np.sqrt(eph['RA_3sigma'][0]**2+eph['DEC_3sigma'][0]**2)

            if ztf_entry_unc < max_unc:
                print('hit: object {}'.format(obj_name))
                
                print('object name:', obj_name)
                print('image date (UT):', ztf_entry_dateUT)
                print('fracday:', ztf_entry_fracday)
                print('filter code:', ztf_entry_filtercode)
                print('distance to the Sun:', ztf_entry_sundist, 'AU')
                print('distance to the Earth:', ztf_entry_geodist, 'AU')
                print('predicted V magnitude:', ztf_entry_vmag)
                print('ecliptic latitude:', eph['ObsEclLat'][0])
                print('positional uncertainty: %.2f arcsec' % ztf_entry_unc)
                
                break
hit: object 2018 CL
object name: 2018 CL
image date (UT): 20180205
fracday: 231528
filter code: zr
distance to the Sun: 0.9936 AU
distance to the Earth: 0.0077 AU
predicted V magnitude: 15.69
ecliptic latitude: -7.4574757
positional uncertainty: 2.31 arcsec

Voila! You made it. Note that it takes quite a bit of time (a day or two) to loop over the entire NEA catalog -- this is because Horizons has to calculate the ephemerides for each of these ~15,000 NEAs and this is slow! Therefore, we provide a pre-generated result file for you to proceed to step 2.

But do take a few known asteroids and see what you get for them.

Step 2: compare the catalog derived from step 1 to the source catalog generated by ZTF streak detection pipeline

Let's read in the result of step 1... (OK, we cheated, it is pre-generated)

In [14]:
check_fmo = np.genfromtxt('check_fmo.txt', dtype=None, \
                          delimiter=(12, 10, 6, 10, 10, 10, 5, 10, 10, 10), \
                         names=["object", "dateUT", "filter", "rate", "sundist", \
                                "geodist", "vmag", "ecllat", "unc", "fracday"])
/home/chummels/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default.
  after removing the cwd from sys.path.

And also we read the source files from ZTF... essentially, these are a large number of "streak_qa" files, each file contains candidate streaks in one image. Apparently, we want to loop over them and see if there is any candidate that coincides the position of a known streaking NEA. We also record the rate and mag of detections and non-detections for diagnostic purpose.

For demonstrative purpose I am picking one of them to show you how this whole thing works. Let us start with our old friend, 2018 CL.

In [15]:
for cfi in check_fmo:
    if cfi[0].strip().decode('UTF8') == '2018 CL':
        break

What is the uncertainty range for this entry?

In [16]:
print(cfi[8], "arcsec")
2.3094 arcsec

This is very precise! Well, we know it should be small because the orbit from JPL has included ZTF observations.

Now, let's obtain the predicted coordinate for 2018 CL at this time, using the knowledge we just learned in Step 1.

In [17]:
expTimeJD = Time('{}-{}-{}'.format(str(cfi[1])[0:4], str(cfi[1])[4:6], str(cfi[1])[6:8]), format='iso').jd + float(cfi[9])/1000000

expTimeJD
Out[17]:
2458154.731528
In [18]:
obj_this = Horizons(id=obj_name, location='I41', epochs=expTimeJD)
eph = obj_this.ephemerides()

eph
Out[18]:
Table masked=True length=1
targetnamedatetime_strdatetime_jdHGsolar_presenceflagsRADECRA_appDEC_appRA_rateDEC_rateAZELAZ_rateEL_ratesat_Xsat_Ysat_PANGsiderealtimeairmassmagextinctVilluminationillum_defectsat_sepsat_visang_widthPDObsLonPDObsLatPDSunLonPDSunLatSubSol_angSubSol_distNPole_angNPole_distEclLonEclLatrr_ratedeltadelta_ratelighttimevel_sunvel_obselongelongFlagalphalunar_elonglunar_illumsat_alphasunTargetPAvelocityPAOrbPlaneAngconstellationTDB-UTObsEclLonObsEclLatNPole_RANPole_DECGlxLonGlxLatsolartimeearth_lighttimeRA_3sigmaDEC_3sigmaSMAA_3sigmaSMIA_3sigmaTheta_3sigmaArea_3sigmaRSS_3sigmar_3sigmar_rate_3sigmaSBand_3sigmaXBand_3sigmaDoppDelay_3sigmatrue_anomhour_anglealpha_truePABLonPABLat
------dmag---------degdegdegdegarcsec / harcsec / hdegdegarcsec / minarcsec / minarcsecarcsecdeg------magmag%arcsecarcsec---arcsecdegdegdegdegdegarcsecdegarcsecdegdegAUkm / sAUkm / sminkm / skm / sdeg---degdeg%degdegdegdeg---sdegdegdegdegdegdeg---minarcsecarcsecarcsecarcsecdegarcsec2arcseckmkm / sHzHzsdeg---degdegdeg
str9str24float64float64float64str1str1float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64int64float64str1int64int64int64int64int64float64float64int64int64float64float64float64float64float64float64float64float64float64float64str2float64float64float64float64float64float64float64str3float64float64float64int64int64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64
(2018 CL)2018-Feb-05 05:33:24.0192458154.73152799825.50.15130.2430410.59115130.4917110.52376-2179.933843.943124.290355.1746612.91693.91594314.595242.96122.7336.78755582461.2170.2115.5599.291--612957.6*----------235.370.0----136.0615-0.05880.993552869965-4.89240910.00768971056396-5.55238570.06395327.287358.81488170.266/T9.658769.673.20.071455.334274.9696.00716Cnc69.18491130.0804119-7.4577146----215.72888129.09429621.53361001240.0003541.321.8952.0980.965-61.0812.722.3091487.72610.0070343111.58405.530.009925217.3296-1.9118911799.6624132.954-3.7651

A simplistic approach is to loop through each small folder generated by findStreaks and see if there is any candidate matching the predicted RA and DEC:

In [19]:
detectionRadius = 5/3600
In [20]:
d = sorted(glob.glob('{}/ztf*'.format('run_merged')))
targetRA = eph['RA'][0]
targetDec = eph['DEC'][0]

for di in d:
    fileName = os.path.basename(di)
    streaksQA_path = '{}/{}_streaksqa.txt'.format(di, fileName)
    
    if os.stat(streaksQA_path).st_size > 500:
        try:
            streaksQA = np.loadtxt(streaksQA_path, delimiter=',', comments='\\')
            for streaksQA_item in streaksQA:
                if (abs(streaksQA_item[1]-targetRA) < detectionRadius and \
                abs(streaksQA_item[2]-targetDec) < detectionRadius) or \
                (abs(streaksQA_item[3]-targetRA) < detectionRadius and \
                abs(streaksQA_item[4]-targetDec) < detectionRadius):
                        print('YES! Streak found in', streaksQA_path)
        except:
            continue
YES! Streak found in run_merged/ztf_20180205231528_000517_zr_c06_o_q4/ztf_20180205231528_000517_zr_c06_o_q4_streaksqa.txt

Wonderful! We have one hit! We can even take a look at the cutout to see what it looks like:

In [21]:
for di in d:
    fileName = os.path.basename(di)
    
    streaksQA_path = '{}/{}_streaksqa.txt'.format(di, fileName)
    
    if os.stat(streaksQA_path).st_size > 500:
        try:
            streaksQA = np.loadtxt(streaksQA_path, delimiter=',', comments='\\')
            for streaksQA_item in streaksQA:
                if (abs(streaksQA_item[1]-targetRA) < detectionRadius and \
                    abs(streaksQA_item[2]-targetDec) < detectionRadius) or \
                    (abs(streaksQA_item[3]-targetRA) < detectionRadius and \
                    abs(streaksQA_item[4]-targetDec) < detectionRadius):
                        folder_name = di.split('/')[-1]
                        jpgs = glob.glob('run_merged/' + folder_name + '/' + folder_name + '_cutouts/*.jpg')
                        display(Image(filename = jpgs[0], width=157, height=157))
        except:
            continue

Voila! Now, homework for you: can you use the example shown above to work out the detection efficiency of ZTF?