search

B612-Asteroid-Institute / precovery

Fast precovery of small body observations at scale
BSD 3-Clause "New" or "Revised" License
5 stars 2 forks source link

Precovery of (458508) #4

Closed moeyensj closed 2 years ago

moeyensj commented 2 years ago

Brute-force attribution of (458508) finds 10 observations within 1 arcsecond in observations dated between August 30th, 2013 and September 30th, 2013.

The precovery code finds 2 observations, only observation c4d.230825.30.319 agrees with the ones found by brute-force attribution.

import os
import pandas as pd 
from astropy.time import Time

from precovery import precovery_db
from precovery.orbit import Orbit
from precovery.orbit import EpochTimescale

SYSTEM_DIR = "/epyc"
ORBIT_DIR = os.path.join(SYSTEM_DIR, "projects/thor/thor_results/nsc/precovery/458508/")
DB_DIR = os.path.join(SYSTEM_DIR, "projects/thor/thor_data/nsc/precovery/")
OBSERVATIONS_FILE =  os.path.join(SYSTEM_DIR, "projects/thor/thor_data/nsc/preprocessed/nsc_dr2_observations_2013-08-30_2013-09-30.h5")

# Read database from directory
db = precovery_db.PrecoveryDatabase.from_dir(DB_DIR, create=True)

# Read MPC (truth) orbit
mpc_orbit = pd.read_csv(
    os.path.join(ORBIT_DIR, "mpc_orbit_cartesian.csv"),
    index_col=False
)
# Read results from brute-force attributions
orbit_attributions = pd.read_csv(
    os.path.join(ORBIT_DIR, "attributions.csv"),
    index_col=False
)
# Read observations (not necessary to reproduce example). 
# observations = pd.read_hdf(OBSERVATIONS_FILE, key="data")

# Create precovery Orbit
# Epochs are defined in TDB, observations are defined in UTC. So lets convert the orbit's 
# epoch to UTC so that when propagation and ephemeris generation occurs, the times scales match
# the observations.
epoch_tdb = Time(mpc_orbit["mjd_tdb"].values[0], scale="tdb", format="mjd")
epoch_utc = epoch_tdb.utc.mjd
orbit = Orbit.cartesian(
    0,
    mpc_orbit["x"].values[0],
    mpc_orbit["y"].values[0],
    mpc_orbit["z"].values[0],
    mpc_orbit["vx"].values[0],
    mpc_orbit["vy"].values[0],
    mpc_orbit["vz"].values[0],
    epoch_utc,
    EpochTimescale.UTC,
    20,
    0.15
)

matches = [m for m in db.precover(orbit, tolerance=1/3600)]
obs_ids_precovered = [m.id for m in matches]
print(obs_ids_precovered)

print(orbit_attributions["obs_id"].values)

This outputs the following:

['c4d.230825.30.319', 'c4d.237746.22.451']
['c4d.230114.59.331' 'c4d.230825.30.319' 'c4d.231493.14.573'
 'c4d.232341.24.965' 'c4d.232342.24.1141' 'c4d.232755.5.890'
 'c4d.232756.5.1516' 'c4d.232759.5.1678' 'c4d.233134.54.352'
 'c4d.233473.32.288']
spenczar commented 2 years ago

Hm, yeah, I'm able to reproduce. Here's some debug-level logging: output.log

spenczar commented 2 years ago

Brute force finds observations at these epochs:

[56538.20262244 56540.05018224 56543.12915145 56545.18994486
 56545.19133431 56546.16615175 56546.16752796 56546.17171658
 56547.17895786 56548.11889986

Each of those should correspond to one or more Frames in the precovery DB. Let's check the log:

56538.20262244:

The window midpoint was <Ephemeris ra=325.2442 dec=-1.1821 mjd=56541.000114> (healpixel: 4751), which is 3 days away.

DEBUG:precovery:mjd=56538.202622:   healpixels with data: {4772, 4773, 4750, 4751, 4784}
DEBUG:precovery:mjd=56538.202622:   ephemeris at ra=328.547 dec=-0.168  healpix=4762

5 pixels with data; we don't hit any of them.

56540.05018224

This is one of the ones we successfully find.

56543.12915145

Midpoint: <Ephemeris ra=325.2442 dec=-1.1821 mjd=56541.000114> (healpixel: 4751)

DEBUG:precovery:mjd=56543.129151:   healpixels with data: {4772, 4773, 4750, 4751}
DEBUG:precovery:mjd=56543.129151:   ephemeris at ra=322.729 dec=-0.222  healpix=4774

Off by one healpixel!

56545.18994486

56545.19133431

We get pretty close on these two; we get a healpix rough match just a bit off:

DEBUG:precovery:mjd=56545.194126:       healpixels with data: {4772, 4773, 4750, 4751, 4784}                                                                                                                                                                                                                                                                
DEBUG:precovery:mjd=56545.194126:       ephemeris at ra=327.502 dec=-0.139      healpix=4762                                                                                                                                                                                                                                                                
DEBUG:precovery:mjd=56545.192749:       healpixels with data: {4760, 4762, 4749, 4751}                                                                                                                                                                                                                                                                      
DEBUG:precovery:mjd=56545.192749:       ephemeris at ra=327.504 dec=-0.139      healpix=4762                                                                                                                                                                                                                                                                
DEBUG:precovery:mjd=56545.192749: healpixel collision, checking frames                                                                                                                                                                                                                                                                                      
INFO:precovery:checking frames for healpix=4762 obscode=W84 mjd=56545.192749                                                                                                                                                                                                                                                                                
INFO:precovery:checking frame: HealpixFrame(id=5708, obscode='W84', catalog_id='c4d_130910_043733_ooi_r_ls9', mjd=56545.19274945464, healpixel=4762, data_uri='frames_00000000.data', data_offset=366546274, data_length=7337)                                                                                                                              
INFO:precovery:checked 199 observations in frame                                                                                                                                                                                                                                                                                                            
INFO:precovery:checked 1 frames                                                                                                                                                                                                                                                                                                                             
DEBUG:precovery:mjd=56545.191334:       healpixels with data: {4769, 4772, 4773, 4750}                                                                                                                                                                                                                                                                      
DEBUG:precovery:mjd=56545.191334:       ephemeris at ra=327.506 dec=-0.139      healpix=4762                                                                                                                                                                                                                                                                
DEBUG:precovery:mjd=56545.189945:       healpixels with data: {4769, 4772, 4773, 4750}                                                                                                                                                                                                                                                                      
DEBUG:precovery:mjd=56545.189945:       ephemeris at ra=327.508 dec=-0.139      healpix=4762                                                                                                                                                                                                                                                                
DEBUG:precovery:mjd=56545.187351:       healpixels with data: {11652, 11653, 11566, 11567}                                                                                                                                                                                                                                                                  
DEBUG:precovery:mjd=56545.187351:       ephemeris at ra=327.511 dec=-0.139      healpix=4762

56546.16615175

56546.16752796

56546.17171658

We barely miss all three of these - we check the actual frames, but don't think it matches:

DEBUG:precovery:mjd=56546.171717:       healpixels with data: {4772, 4773, 4747, 4750, 4751}                       
DEBUG:precovery:mjd=56546.171717:       ephemeris at ra=326.265 dec=-0.152      healpix=4751                       
DEBUG:precovery:mjd=56546.171717: healpixel collision, checking frames                                             
INFO:precovery:checking frames for healpix=4751 obscode=W84 mjd=56546.171717                                       
INFO:precovery:checking frame: HealpixFrame(id=6256, obscode='W84', catalog_id='c4d_130911_040716_ooi_i_d2', mjd=5\
6546.17171658296, healpixel=4751, data_uri='frames_00000000.data', data_offset=385680752, data_length=1794)        
INFO:precovery:checked 49 observations in frame                                                                    
INFO:precovery:checked 1 frames                                                                                    
DEBUG:precovery:mjd=56546.170297:       healpixels with data: {4768, 4769, 4770, 4771}                             
DEBUG:precovery:mjd=56546.170297:       ephemeris at ra=326.267 dec=-0.152      healpix=4751                       
DEBUG:precovery:mjd=56546.168919:       healpixels with data: {4768, 4769, 4770, 4771}                             
DEBUG:precovery:mjd=56546.168919:       ephemeris at ra=326.269 dec=-0.152      healpix=4751                       
DEBUG:precovery:mjd=56546.167528:       healpixels with data: {4772, 4773, 4747, 4750, 4751}                       
DEBUG:precovery:mjd=56546.167528:       ephemeris at ra=326.271 dec=-0.152      healpix=4751                       
DEBUG:precovery:mjd=56546.167528: healpixel collision, checking frames                                             
INFO:precovery:checking frames for healpix=4751 obscode=W84 mjd=56546.167528                                       
INFO:precovery:checking frame: HealpixFrame(id=6243, obscode='W84', catalog_id='c4d_130911_040114_ooi_r_ls9', mjd=\
56546.16752795596, healpixel=4751, data_uri='frames_00000000.data', data_offset=385120817, data_length=2981)       
INFO:precovery:checked 81 observations in frame                                                                    
INFO:precovery:checked 1 frames                                                                                    
DEBUG:precovery:mjd=56546.166152:       healpixels with data: {4772, 4773, 4750, 4751}                             
DEBUG:precovery:mjd=56546.166152:       ephemeris at ra=326.272 dec=-0.152      healpix=4751                       
DEBUG:precovery:mjd=56546.166152: healpixel collision, checking frames                                             
INFO:precovery:checking frames for healpix=4751 obscode=W84 mjd=56546.166152                                       
INFO:precovery:checking frame: HealpixFrame(id=6238, obscode='W84', catalog_id='c4d_130911_035915_ooi_g_ls9', mjd=\
56546.16615174757, healpixel=4751, data_uri='frames_00000000.data', data_offset=384942080, data_length=2808)       
INFO:precovery:checked 77 observations in frame                                                                    
INFO:precovery:checked 1 frames

56547.17895786

DEBUG:precovery:mjd=56547.178958:       healpixels with data: {4772, 4773, 4774, 4775}                             
DEBUG:precovery:mjd=56547.178958:       ephemeris at ra=324.988 dec=-0.166      healpix=4773                       
DEBUG:precovery:mjd=56547.178958: healpixel collision, checking frames                                             
INFO:precovery:checking frames for healpix=4773 obscode=W84 mjd=56547.178958                                       
INFO:precovery:checking frame: HealpixFrame(id=7311, obscode='W84', catalog_id='c4d_130912_041741_ooi_g_ls9', mjd=\
56547.17895786278, healpixel=4773, data_uri='frames_00000000.data', data_offset=403607979, data_length=22002)      
INFO:precovery:checked 601 observations in frame                                                                   
INFO:precovery:checked 1 frames

56548.11889986

DEBUG:precovery:mjd=56548.118900:       healpixels with data: {4772, 4773, 4750, 4751}                             
DEBUG:precovery:mjd=56548.118900:       ephemeris at ra=323.796 dec=-0.179      healpix=4773                       
DEBUG:precovery:mjd=56548.118900: healpixel collision, checking frames                                             
INFO:precovery:checking frames for healpix=4773 obscode=W84 mjd=56548.118900                                       
INFO:precovery:checking frame: HealpixFrame(id=8041, obscode='W84', catalog_id='c4d_130913_025112_ooi_i_d2', mjd=5\
6548.11889985669, healpixel=4773, data_uri='frames_00000000.data', data_offset=426112732, data_length=45828)       
INFO:precovery:checked 1236 observations in frame                                                                  
INFO:precovery:checked 1 frames

Summary

Out of 10:

spenczar commented 2 years ago

Turned on more extensive logging, which lets me see the ID of every observation that the precovery scan touched. It didn't examine any of the ones that brute force found.

This means we are definitely skipping them based on healpix calculations.

spenczar commented 2 years ago

@moeyensj Can you provide the RA and Dec of the brute force ephemeris, so we can compare them to the ones that show up in these logs?

I wonder if the approximate propagation is too inaccurate, placing us in the wrong healpixel.

We may need to use coarser healpixels, if so.

moeyensj commented 2 years ago

Looping back here: I just ran the latest version of precovery with this particular case and we now recover the same observations as the brute-force attribution code (and more with the full dataset). 

orbit_id | mjd_utc | ra_deg | dec_deg | ra_sigma_arcsec | dec_sigma_arcsec | mag | mag_sigma | filter | obscode | exposure_id | observation_id | healpix_id | pred_ra_deg | pred_dec_deg | pred_vra_degpday | pred_vdec_degpday | delta_ra_arcsec | delta_dec_arcsec | distance_arcsec
458508 (2011 CM42) | 56538.202622 | 325.802459 | -1.153522 | 0.026354 | 0.025239 | 20.947607 | 0.028126 | g | W84 | c4d_130903_045146_ooi_g_ls9 | c4d.230114.59.331 | 4982414964 | 325.802563 | -1.153512 | -0.208044 | -0.009353 | 0.374038 | 0.036238 | 0.375714
458508 (2011 CM42) | 56540.050182 | 325.431560 | -1.171906 | 0.033156 | 0.032910 | 21.119007 | 0.030525 | g | W84 | c4d_130905_011215_ooi_g_ls9 | c4d.230825.30.319 | 4982463327 | 325.431671 | -1.171895 | -0.202804 | -0.010531 | 0.398636 | 0.040449 | 0.400600
458508 (2011 CM42) | 56543.129151 | 324.832448 | -1.206604 | 0.044687 | 0.051652 | 20.568155 | 0.040078 | g | W84 | c4d_130908_030558_ooi_g_ls9 | c4d.231493.14.573 | 4981784508 | 324.832553 | -1.206565 | -0.194785 | -0.011956 | 0.380399 | 0.139689 | 0.405157
458508 (2011 CM42) | 56545.189945 | 324.447425 | -1.232025 | 0.026480 | 0.027871 | 20.908384 | 0.036654 | r | W84 | c4d_130910_043331_ooi_r_ls9 | c4d.232341.24.965 | 5004147711 | 324.447540 | -1.232011 | -0.187796 | -0.012658 | 0.414371 | 0.048976 | 0.417160
458508 (2011 CM42) | 56545.191334 | 324.447176 | -1.232045 | 0.027617 | 0.030041 | 20.741930 | 0.038616 | i | W84 | c4d_130910_043531_ooi_i_d2 | c4d.232342.24.1141 | 5004147711 | 324.447280 | -1.232029 | -0.187775 | -0.012658 | 0.371052 | 0.057621 | 0.375414
458508 (2011 CM42) | 56546.166152 | 324.270161 | -1.244569 | 0.028601 | 0.035915 | 20.934050 | 0.026578 | g | W84 | c4d_130911_035915_ooi_g_ls9 | c4d.232755.5.890 | 5004168139 | 324.270277 | -1.244556 | -0.184676 | -0.012954 | 0.414544 | 0.045944 | 0.416985
458508 (2011 CM42) | 56546.167528 | 324.269930 | -1.244591 | 0.023201 | 0.024901 | 20.651749 | 0.023146 | r | W84 | c4d_130911_040114_ooi_r_ls9 | c4d.232756.5.1516 | 5004168139 | 324.270022 | -1.244574 | -0.184661 | -0.012954 | 0.333352 | 0.063443 | 0.339258
458508 (2011 CM42) | 56546.171717 | 324.269153 | -1.244645 | 0.025260 | 0.027846 | 20.624908 | 0.031053 | i | W84 | c4d_130911_040716_ooi_i_d2 | c4d.232759.5.1678 | 5004168139 | 324.269249 | -1.244628 | -0.184611 | -0.012952 | 0.344765 | 0.063469 | 0.350478
458508 (2011 CM42) | 56547.178958 | 324.089568 | -1.257856 | 0.067939 | 0.067976 | 21.521446 | 0.055080 | g | W84 | c4d_130912_041741_ooi_g_ls9 | c4d.233134.54.352 | 5004131115 | 324.089666 | -1.257840 | -0.180970 | -0.013194 | 0.352896 | 0.056755 | 0.357346
458508 (2011 CM42) | 56548.118900 | 323.925554 | -1.270400 | 0.026000 | 0.028242 | 20.773983 | 0.033952 | i | W84 | c4d_130913_025112_ooi_i_d2 | c4d.233473.32.288 | 5004216738 | 323.925645 | -1.270377 | -0.177881 | -0.013421 | 0.326448 | 0.082351 | 0.336598
458508 (2011 CM42) | 56548.121702 | 323.925048 | -1.270425 | 0.047149 | 0.047211 | 20.465689 | 0.058700 | z | W84 | c4d_130913_025515_ooi_z_ls9 | c4d.233475.32.268 | 5004216738 | 323.925146 | -1.270415 | -0.177870 | -0.013419 | 0.354843 | 0.034435 | 0.356424
458508 (2011 CM42) | 56558.175074 | 322.383770 | -1.407888 | 0.075327 | 0.075890 | 21.017221 | 0.074004 | z | W84 | c4d_130923_041206_ooi_z_ls9 | c4d.237746.22.451 | 5004463130 | 322.383846 | -1.407892 | -0.136162 | -0.013201 | 0.272960 | -0.014053 | 0.273239
458508 (2011 CM42) | 56559.040490 | 322.271682 | -1.419309 | 0.029465 | 0.032234 | 20.957930 | 0.027108 | g | W84 | c4d_130924_005818_ooi_g_ls9 | c4d.237968.45.115 | 5001671998 | 322.271670 | -1.419296 | -0.132839 | -0.013095 | -0.044187 | 0.047858 | 0.065128
458508 (2011 CM42) | 56562.105627 | 321.898999 | -1.458090 | 0.074000 | 0.073323 | 21.831330 | 0.067864 | g | W84 | c4d_130927_023206_ooi_g_a1 | c4d.238654.14.643 | 5001681916 | 321.899054 | -1.458074 | -0.118946 | -0.012137 | 0.196576 | 0.059472 | 0.205315
458508 (2011 CM42) | 56562.106660 | 321.898713 | -1.458106 | 0.032952 | 0.037530 | 21.746834 | 0.027630 | g | W84 | c4d_130927_023335_ooi_g_a1 | c4d.238655.20.1338 | 5001681916 | 321.898931 | -1.458086 | -0.118936 | -0.012136 | 0.785701 | 0.070813 | 0.788632
458508 (2011 CM42) | 56562.110465 | 321.898269 | -1.458142 | 0.028482 | 0.034508 | 21.859716 | 0.024440 | g | W84 | c4d_130927_023904_ooi_g_a1 | c4d.238656.20.1390 | 5001681916 | 321.898479 | -1.458133 | -0.118894 | -0.012132 | 0.754359 | 0.033398 | 0.754854
458508 (2011 CM42) | 56562.121918 | 321.897063 | -1.458258 | 0.027734 | 0.028268 | 20.976082 | 0.020990 | i | W84 | tu1870287 | c4d.238659.20.2113 | 5001681915 | 321.897118 | -1.458272 | -0.118755 | -0.012120 | 0.195874 | -0.050717 | 0.202272
458508 (2011 CM42) | 56562.125720 | 321.896474 | -1.458344 | 0.040462 | 0.040975 | 21.497398 | 0.034360 | i | W84 | tu1869880 | c4d.238660.20.2023 | 5001681915 | 321.896666 | -1.458318 | -0.118704 | -0.012116 | 0.692527 | 0.095736 | 0.698890
458508 (2011 CM42) | 56575.077207 | 320.815386 | -1.573892 | 0.021589 | 0.022694 | 21.144255 | 0.032908 | r | W84 | c4d_131010_015110_ooi_r_ls9 | c4d.242417.45.769 | 5002970670 | 320.815418 | -1.573873 | -0.056226 | -0.004746 | 0.114876 | 0.069282 | 0.134114
458508 (2011 CM42) | 56577.014406 | 320.724020 | -1.581685 | 0.066623 | 0.066406 | 21.874924 | 0.124678 | z | W84 | c4d_131012_002044_ooi_z_ls9 | c4d.243034.14.1925 | 5002960755 | 320.724033 | -1.581687 | -0.046872 | -0.003265 | 0.047167 | -0.005262 | 0.047441
458508 (2011 CM42) | 57827.325168 | 168.184770 | -30.171123 | 0.023761 | 0.025492 | 19.466130 | 0.020774 | i | W84 | c4d_170315_075109_ooi_i_v1 | c4d.630397.5.1748 | 10155775404 | 168.184978 | -30.171095 | -0.302032 | -0.050202 | 0.746321 | 0.100066 | 0.652930
458508 (2011 CM42) | 57830.257999 | 167.320653 | -30.303282 | 0.017331 | 0.014685 | 19.291792 | 0.013374 | i | W84 | c4d_170318_061352_ooi_i_v1 | c4d.631391.57.555 | 10155704895 | 167.320801 | -30.303279 | -0.300780 | -0.034751 | 0.535068 | 0.008931 | 0.462046