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
Compute the detection rate of fast-moving asteroids, asteroids that pass within ~20 lunar distances from Earth, and leave trails/streaks on astronomical images.
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
).
None
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
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).
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:
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',)
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.
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...
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?
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:
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.
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?
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.
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".
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!)
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.
Let's read in the result of step 1... (OK, we cheated, it is pre-generated)
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.
for cfi in check_fmo:
if cfi[0].strip().decode('UTF8') == '2018 CL':
break
What is the uncertainty range for this entry?
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.
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
2458154.731528
obj_this = Horizons(id=obj_name, location='I41', epochs=expTimeJD)
eph = obj_this.ephemerides()
eph
targetname | datetime_str | datetime_jd | H | G | solar_presence | flags | RA | DEC | RA_app | DEC_app | RA_rate | DEC_rate | AZ | EL | AZ_rate | EL_rate | sat_X | sat_Y | sat_PANG | siderealtime | airmass | magextinct | V | illumination | illum_defect | sat_sep | sat_vis | ang_width | PDObsLon | PDObsLat | PDSunLon | PDSunLat | SubSol_ang | SubSol_dist | NPole_ang | NPole_dist | EclLon | EclLat | r | r_rate | delta | delta_rate | lighttime | vel_sun | vel_obs | elong | elongFlag | alpha | lunar_elong | lunar_illum | sat_alpha | sunTargetPA | velocityPA | OrbPlaneAng | constellation | TDB-UT | ObsEclLon | ObsEclLat | NPole_RA | NPole_DEC | GlxLon | GlxLat | solartime | earth_lighttime | RA_3sigma | DEC_3sigma | SMAA_3sigma | SMIA_3sigma | Theta_3sigma | Area_3sigma | RSS_3sigma | r_3sigma | r_rate_3sigma | SBand_3sigma | XBand_3sigma | DoppDelay_3sigma | true_anom | hour_angle | alpha_true | PABLon | PABLat |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
--- | --- | d | mag | --- | --- | --- | deg | deg | deg | deg | arcsec / h | arcsec / h | deg | deg | arcsec / min | arcsec / min | arcsec | arcsec | deg | --- | --- | mag | mag | % | arcsec | arcsec | --- | arcsec | deg | deg | deg | deg | deg | arcsec | deg | arcsec | deg | deg | AU | km / s | AU | km / s | min | km / s | km / s | deg | --- | deg | deg | % | deg | deg | deg | deg | --- | s | deg | deg | deg | deg | deg | deg | --- | min | arcsec | arcsec | arcsec | arcsec | deg | arcsec2 | arcsec | km | km / s | Hz | Hz | s | deg | --- | deg | deg | deg |
str9 | str24 | float64 | float64 | float64 | str1 | str1 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | int64 | float64 | str1 | int64 | int64 | int64 | int64 | int64 | float64 | float64 | int64 | int64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | str2 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | str3 | float64 | float64 | float64 | int64 | int64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
(2018 CL) | 2018-Feb-05 05:33:24.019 | 2458154.731527998 | 25.5 | 0.15 | 130.24304 | 10.59115 | 130.49171 | 10.52376 | -2179.93 | 3843.943 | 124.2903 | 55.1746 | 612.91 | 693.91 | 594314.5 | 95242.96 | 122.733 | 6.7875558246 | 1.217 | 0.21 | 15.55 | 99.291 | -- | 612957.6 | * | -- | -- | -- | -- | -- | 235.37 | 0.0 | -- | -- | 136.0615 | -0.0588 | 0.993552869965 | -4.8924091 | 0.00768971056396 | -5.5523857 | 0.063953 | 27.28735 | 8.81488 | 170.266 | /T | 9.6587 | 69.6 | 73.2 | 0.0714 | 55.334 | 274.969 | 6.00716 | Cnc | 69.18491 | 130.0804119 | -7.4577146 | -- | -- | 215.728881 | 29.094296 | 21.5336100124 | 0.000354 | 1.32 | 1.895 | 2.098 | 0.965 | -61.08 | 12.72 | 2.309 | 1487.7261 | 0.0070343 | 111.58 | 405.53 | 0.009925 | 217.3296 | -1.911891179 | 9.6624 | 132.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:
detectionRadius = 5/3600
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:
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?