Merge pull request #97 from rahil-makadia/dev

Minor updates

Author: rahil-makadia

grss/debias/get_debiasing_data.py

@@ -8,12 +8,12 @@
 # if lowres_data.tgz or lowres_data/ already exist, do not download
 if not os.path.exists(f'{cwd}/lowres_data.tgz') and not os.path.exists(f'{cwd}/lowres_data/'):
     FTP_FILE = 'debias_2018.tgz'
-    cmd = f'curl --connect-timeout 30 --silent --show-error -o {cwd}/lowres_data.tgz {FTP_SITE}/{FTP_FILE}'
+    cmd = f'curl --connect-timeout 30 --show-error -o {cwd}/lowres_data.tgz {FTP_SITE}/{FTP_FILE}'
     print('Downloading low-resolution debiasing data...')
     os.system(cmd)
 if not os.path.exists(f'{cwd}/hires_data.tgz') and not os.path.exists(f'{cwd}/hires_data/'):
     FTP_FILE = 'debias_hires2018.tgz'
-    cmd = f'curl --connect-timeout 30 --silent --show-error -o {cwd}/hires_data.tgz {FTP_SITE}/{FTP_FILE}'
+    cmd = f'curl --connect-timeout 30 --show-error -o {cwd}/hires_data.tgz {FTP_SITE}/{FTP_FILE}'
     print('Downloading high-resolution debiasing data...')
     os.system(cmd)
 

grss/fit/fit_ades.py

@@ -110,8 +110,8 @@
 
 special_codes = {
     'gaia': {'258'},
-    'occultation': {'275'},
-    'spacecraft': {'S/C', 'S_C', '245', '249', '250', '274', 'C49', 'C50', 'C51',
-                   'C52', 'C53', 'C54', 'C55', 'C56', 'C57', 'C59', },
+    'occultation': {'244', '275'},
+    'spacecraft': {'S/C', 'S_C', '245', '249', '250', '273', '274',
+                   'C49', 'C50', 'C51', 'C52', 'C53', 'C54', 'C55', 'C56', 'C57', 'C58', 'C59', },
     'roving': {'247', '270'},
 }

grss/fit/fit_optical.py

@@ -67,8 +67,6 @@ def _ades_mode_check(df):
     if not df['mode'].isin(valid_modes).all():
         raise ValueError(f"Invalid mode in the data: {df['mode'].unique()}.\n"
                         f"Acceptable modes are {valid_modes}.")
-    if df['mode'].isin(['UNK']).any():
-        print("\tWARNING: At least one unknown observation mode in the data.")
     return None
 
 def _ades_ast_cat_check(df):
@@ -102,8 +100,6 @@ def _ades_ast_cat_check(df):
         'ACRS', 'LickGas', 'Ida93', 'Perth70', 'COSMOS',
         'Yale', 'ZZCAT', 'IHW', 'GZ', 'UNK']
     df_cats = df['astCat']
-    if df_cats.isin(deprecated_cats).any():
-        print("\tWARNING: At least one deprecated star catalog in the data.")
     if not df_cats.isin(valid_cats).all():
         invalid_cats = np.setdiff1d(df_cats.unique(), valid_cats)
         print("\tWARNING: At least one unrecognized astCat in the data. "
@@ -218,11 +214,24 @@ def add_psv_obs(psv_obs_file, obs_df, t_min_tdb=None, t_max_tdb=None, verbose=Fa
         t_min_tdb = -np.inf
     if t_max_tdb is None:
         t_max_tdb = np.inf
-    psv_df = pd.read_csv(psv_obs_file, sep='|', skipinitialspace=True)
+    # read file and skip header rows starting with '#' or '!'
+    skip_count = 0
+    with open(psv_obs_file, 'r') as f:
+        lines = f.readlines()
+        for line in lines:
+            if line.startswith('#') or line.startswith('!'):
+                skip_count += 1
+            else:
+                break
+    psv_df = pd.read_csv(psv_obs_file, sep='|', skipinitialspace=True, skiprows=skip_count)
     psv_df.columns = psv_df.columns.str.strip()
     psv_df = psv_df[[col for col in psv_df.columns if col in ades_column_types]]
     psv_column_types = {col: ades_column_types[col] for col in psv_df.columns}
     psv_df = psv_df.astype(psv_column_types)
+    psv_df['permID'] = obs_df.iloc[-1]['permID']
+    psv_df['provID'] = obs_df.iloc[-1]['provID']
+    # strip every entry in a column of leading and trailing whitespace
+    psv_df = psv_df.apply(lambda x: x.str.strip() if x.dtype == "object" else x)
     occ_idx = psv_df.query("mode == 'OCC'").index
     if len(occ_idx) > 0:
         psv_df.loc[occ_idx,'dec'] = psv_df.loc[occ_idx,'decStar']
@@ -238,7 +247,10 @@ def add_psv_obs(psv_obs_file, obs_df, t_min_tdb=None, t_max_tdb=None, verbose=Fa
     psv_df['cosDec'] = np.cos(psv_df['dec']*np.pi/180)
     psv_df['sigRA'] = psv_df['rmsRA']
     psv_df['sigDec'] = psv_df['rmsDec']
-    psv_df['sigCorr'] = psv_df['rmsCorr']
+    if 'rmsCorr' not in psv_df:
+        psv_df['sigCorr'] = 0.0
+    if 'rmsTime' not in psv_df:
+        psv_df['sigTime'] = 1.0
     times = Time(psv_df['obsTime'].to_list(), format='isot', scale='utc')
     psv_df['obsTimeMJD'] = times.utc.mjd
     psv_df['obsTimeMJDTDB'] = times.tdb.mjd
@@ -287,7 +299,8 @@ def _get_gaia_query_results(body_id, release):
         + "ra_error_systematic,dec_error_systematic,ra_dec_correlation_systematic,"
         + "ra_error_random,dec_error_random,ra_dec_correlation_random,"
         + "x_gaia_geocentric,y_gaia_geocentric,z_gaia_geocentric,"
-        + "vx_gaia_geocentric,vy_gaia_geocentric,vz_gaia_geocentric"
+        + "vx_gaia_geocentric,vy_gaia_geocentric,vz_gaia_geocentric,"
+        + "astrometric_outcome_ccd, astrometric_outcome_transit"
         + f" FROM {release}.{table} {match_str} ORDER BY epoch_utc ASC"
     )
     job = Gaia.launch_job_async(query, dump_to_file=False,background=True)
@@ -298,9 +311,9 @@ def _get_gaia_query_results(body_id, release):
     res.sort('epoch_utc')
     return res
 
-def add_gaia_obs(obs_df, t_min_tdb=None, t_max_tdb=None, gaia_dr='gaiadr3', verbose=False):
+def add_gaia_obs(obs_df, t_min_tdb=None, t_max_tdb=None, gaia_dr='gaiafpr', verbose=False):
     """
-    Assemble the optical observations for a given body from Gaia DR3.
+    Assemble the optical observations for a given body from Gaia FPR.
 
     Parameters
     ----------
@@ -312,7 +325,7 @@ def add_gaia_obs(obs_df, t_min_tdb=None, t_max_tdb=None, gaia_dr='gaiadr3', verb
     t_max_tdb : float, optional
         Maximum time (MJD TDB) for observations to be included, by default None
     gaia_dr : str, optional
-        Gaia data release version database name, by default 'gaiadr3'
+        Gaia data release version database name, by default 'gaiafpr'
     verbose : bool, optional
         Flag to print out information about the observations, by default False
 
@@ -339,6 +352,10 @@ def add_gaia_obs(obs_df, t_min_tdb=None, t_max_tdb=None, gaia_dr='gaiadr3', verb
     curr_transit = int(-1e6)
     gaia_add_counter = 0
     for i, data in enumerate(res):
+        # # print(data["astrometric_outcome_ccd"], data["astrometric_outcome_transit"])
+        # # print(type(data["astrometric_outcome_ccd"]), type(data["astrometric_outcome_transit"]))
+        # if data['astrometric_outcome_ccd'] != 1 or data['astrometric_outcome_transit'] != 1:
+        #     continue
         if curr_transit != data['transit_id']:
             curr_transit = data['transit_id']
             transit_count = 1
@@ -377,16 +394,27 @@ def add_gaia_obs(obs_df, t_min_tdb=None, t_max_tdb=None, gaia_dr='gaiadr3', verb
         obs_df.loc[idx, 'sigRA'] = ra_sig
         obs_df.loc[idx, 'sigDec'] = dec_sig
         obs_df.loc[idx, 'sigCorr'] = corr
+        # obs_df.loc[idx, 'sigTime'] = data['epoch_err']*86400.0
         obs_df.loc[idx, 'biasRA'] = 0.0
         obs_df.loc[idx, 'biasDec'] = 0.0
         obs_df.loc[idx, 'ctr'] = ctr
         obs_df.loc[idx, 'sys'] = sys
-        obs_df.loc[idx, 'pos1'] = data['x_gaia_geocentric']
-        obs_df.loc[idx, 'pos2'] = data['y_gaia_geocentric']
-        obs_df.loc[idx, 'pos3'] = data['z_gaia_geocentric']
-        obs_df.loc[idx, 'vel1'] = data['vx_gaia_geocentric']
-        obs_df.loc[idx, 'vel2'] = data['vy_gaia_geocentric']
-        obs_df.loc[idx, 'vel3'] = data['vz_gaia_geocentric']
+        tcb_tdb_fac = 1 - 1.550519768e-8
+        obs_df.loc[idx, 'pos1'] = data['x_gaia_geocentric']*tcb_tdb_fac
+        obs_df.loc[idx, 'pos2'] = data['y_gaia_geocentric']*tcb_tdb_fac
+        obs_df.loc[idx, 'pos3'] = data['z_gaia_geocentric']*tcb_tdb_fac
+        # obs_df.loc[idx, 'vel1'] = data['vx_gaia_geocentric']*tcb_tdb_fac
+        # obs_df.loc[idx, 'vel2'] = data['vy_gaia_geocentric']*tcb_tdb_fac
+        # obs_df.loc[idx, 'vel3'] = data['vz_gaia_geocentric']*tcb_tdb_fac
+        # # add some position uncertainty (testing)
+        # au2km = 149597870.7
+        # pos_sig = 1.0/au2km
+        # obs_df.loc[idx, 'posCov11'] = pos_sig**2
+        # obs_df.loc[idx, 'posCov12'] = 0.0
+        # obs_df.loc[idx, 'posCov13'] = 0.0
+        # obs_df.loc[idx, 'posCov22'] = pos_sig**2
+        # obs_df.loc[idx, 'posCov23'] = 0.0
+        # obs_df.loc[idx, 'posCov33'] = pos_sig**2
     if verbose:
         print(f"\tFiltered to {gaia_add_counter} observations that",
                 "satisfy the time range constraints.")
@@ -889,10 +917,9 @@ def deweight_obs(obs_df, eff_obs_per_night, verbose):
     times = obs_df['obsTimeMJD'].values
     stations = obs_df['stn'].values
     weights = obs_df[['sigRA', 'sigDec']].values
+    batch_first_time = times[0]
     for i in range(1, len(obs_df)):
-        curr_jd_night = np.floor(times[i]+2400000.5)
-        prev_jd_night = np.floor(times[i-1]+2400000.5)
-        night_match = curr_jd_night == prev_jd_night
+        night_match = times[i] - batch_first_time < 8/24
         curr_observatory = stations[i]
         prev_observatory = stations[i-1]
         observatory_match = curr_observatory == prev_observatory
@@ -904,6 +931,7 @@ def deweight_obs(obs_df, eff_obs_per_night, verbose):
                 factor = night_count**0.5/eff_obs_per_night**0.5
                 weights[i-night_count:i] *= factor
             night_count = 1
+            batch_first_time = times[i]
     # edge case where last observation is part of a batch that needs deweighting
     if night_count > eff_obs_per_night:
         deweight_count += night_count
@@ -939,10 +967,9 @@ def eliminate_obs(obs_df, max_obs_per_night, verbose):
         print(f"Applying {max_obs_per_night}-observation per night elimination scheme.")
     times = obs_df['obsTimeMJD'].values
     stations = obs_df['stn'].values
+    batch_first_time = times[0]
     for i in range(1, len(obs_df)):
-        curr_jd_night = np.floor(times[i]+2400000.5)
-        prev_jd_night = np.floor(times[i-1]+2400000.5)
-        night_match = curr_jd_night == prev_jd_night
+        night_match = times[i] - batch_first_time < 8/24
         curr_observatory = stations[i]
         prev_observatory = stations[i-1]
         observatory_match = curr_observatory == prev_observatory
@@ -953,6 +980,7 @@ def eliminate_obs(obs_df, max_obs_per_night, verbose):
                 idx_to_drop.append(i)
         else:
             night_count = 1
+            batch_first_time = times[i]
     obs_df.drop(idx_to_drop, inplace=True)
     obs_df.reset_index(drop=True, inplace=True)
     if verbose:
@@ -1004,16 +1032,17 @@ def get_optical_obs(body_id, optical_obs_file=None, t_min_tdb=None,
     """
     if eliminate and deweight:
         raise ValueError('Cannot deweight and eliminate observations at the same time.')
-    if not eliminate and not deweight:
+    if not eliminate and not deweight and verbose:
         print("WARNING: No deweighting or elimination scheme applied",
                 "for observations during the same night.")
     obs_df = create_optical_obs_df(body_id, optical_obs_file,
                                     t_min_tdb, t_max_tdb, verbose)
     if debias_lowres is not None:
         obs_df = apply_debiasing_scheme(obs_df, debias_lowres, verbose)
     else:
-        print("WARNING: No debiasing scheme applied to the observations.",
-                "Setting biases to zero.")
+        if verbose:
+            print("WARNING: No debiasing scheme applied to the observations.",
+                    "Setting biases to zero.")
         opt_idx = obs_df.query("mode != 'RAD'").index
         obs_df.loc[opt_idx, ['biasRA', 'biasDec']] = 0.0
     if not accept_weights:

grss/fit/fit_simulation.py

@@ -1366,14 +1366,24 @@ def _inflate_uncertainties(self, prop_sim_past, prop_sim_future):
         sig_ra_vals = self.obs.sigRA.values
         sig_dec_vals = self.obs.sigDec.values
         sig_corr_vals = self.obs.sigCorr.values
+        sys_vals = self.obs.sys.values
+        p11_vals = self.obs.posCov11.values
+        p12_vals = self.obs.posCov12.values
+        p13_vals = self.obs.posCov13.values
+        p22_vals = self.obs.posCov22.values
+        p23_vals = self.obs.posCov23.values
+        p33_vals = self.obs.posCov33.values
         if self.past_obs_exist and self.future_obs_exist:
             optical_obs_dot = prop_sim_past.opticalObsDot + prop_sim_future.opticalObsDot
+            optical_obs_partials = prop_sim_past.opticalPartials + prop_sim_future.opticalPartials
             state_eval = prop_sim_past.xIntegEval + prop_sim_future.xIntegEval
         elif self.past_obs_exist:
             optical_obs_dot = prop_sim_past.opticalObsDot
+            optical_obs_partials = prop_sim_past.opticalPartials
             state_eval = prop_sim_past.xIntegEval
         elif self.future_obs_exist:
             optical_obs_dot = prop_sim_future.opticalObsDot
+            optical_obs_partials = prop_sim_future.opticalPartials
             state_eval = prop_sim_future.xIntegEval
         computed_obs_dot = np.array(optical_obs_dot)[:,0:2]
         computed_obs_dot[:, 0] *= self.obs.cosDec.values
@@ -1399,6 +1409,7 @@ def _inflate_uncertainties(self, prop_sim_past, prop_sim_future):
                 cov = np.array([[sig_ra**2, off_diag],
                                 [off_diag, sig_dec**2]])
                 time_uncert = sig_times[i]
+                # apply time uncertainties to the optical observations
                 if time_uncert != 0.0 and stations[i] not in no_time_uncert:
                     ra_dot_cos_dec = computed_obs_dot[i, 0]
                     dec_dot = computed_obs_dot[i, 1]
@@ -1407,10 +1418,28 @@ def _inflate_uncertainties(self, prop_sim_past, prop_sim_future):
                                         [off_diag_time, dec_dot**2]])*(time_uncert/86400)**2
                     cov += cov_time
                     sig_times[i] = time_uncert
+                # apply Gaia astrometric handling
                 if radius_nonzero and stations[i] == '258':
                     cov_fac = np.array([[fac[i], 0.0],
                                         [0.0, fac[i]]])
                     cov += cov_fac
+                # apply station position uncertainties
+                if np.isfinite(p11_vals[i]):
+                    cov_stn_pos = np.array([
+                        [p11_vals[i], p12_vals[i], p13_vals[i]],
+                        [p12_vals[i], p22_vals[i], p23_vals[i]],
+                        [p13_vals[i], p23_vals[i], p33_vals[i]]])
+                    # make sure all cov values are finite
+                    if not np.isfinite(cov_stn_pos).all():
+                        raise ValueError("Station position covariance matrix is not finite.")
+                    if sys_vals[i] not in {'ICRF_AU', 'ICRF_KM'}:
+                        raise ValueError("Station position covariance matrix system not defined/recognized.")
+                    conv = 1 if sys_vals[i] == 'ICRF_AU' else 1/1.495978707e8
+                    cov_stn_pos *= conv**2
+                    partials = np.array(optical_obs_partials[i]).reshape(2, 6)[:, :3]
+                    partials[0, :] *= self.obs.cosDec.values[i]
+                    cov_stn_unc = partials @ cov_stn_pos @ partials.T
+                    cov += cov_stn_unc
                 det = cov[0, 0]*cov[1, 1] - cov[0, 1]*cov[1, 0]
                 inv = np.array([[cov[1, 1], -cov[0, 1]],
                                 [-cov[1, 0], cov[0, 0]]])/det

grss/fit/fit_utils.py

@@ -92,13 +92,8 @@ def get_mpc_observatory_info():
         dictionary of observatory codes and their corresponding longitude,
         rho*cos(latitude), and rho*sin(latitude)
     """
-    fpath = f'{grss_path}/fit/codes.json'
-    with open(fpath, 'r', encoding='utf-8') as f:
-        codes = json.load(f)
-    mpc_info_dict = {}
-    for key, val in codes.items():
-        if 'Longitude' in val and 'cos' in val and 'sin' in val:
-            mpc_info_dict[key] = val['Longitude'], val['cos'], val['sin']
+    with open(f'{grss_path}/fit/codes.json', 'r', encoding='utf-8') as f:
+        mpc_info_dict = json.load(f)
     return mpc_info_dict
 
 def get_radar_codes_dict():
@@ -287,7 +282,9 @@ def get_sbdb_elems(tdes, cov_elems=True):
         print("get_radar_obs_array: ERROR. ", response.status_code, response.content)
         raise ValueError("Failed to get JPL orbit data")
     raw_data = get_sbdb_raw_data(tdes)
-    if cov_elems:
+    # check that covariance exists
+    cov_exists = raw_data['orbit']['covariance'] is not None
+    if cov_elems and cov_exists:
         # epoch of orbital elements at reference time [JD -> MJD]
         epoch_mjd = float(raw_data['orbit']['covariance']['epoch']) - 2400000.5
         # cometary elements at epoch_mjd
@@ -345,11 +342,16 @@ def get_sbdb_info(tdes):
     # cometary elements corresponding to the covariance on JPL SBDB
     elements = get_sbdb_elems(tdes, cov_elems=True)
     # covariance matrix for cometary orbital elements
-    cov_keys = raw_data['orbit']['covariance']['labels']
-    cov_mat = (np.array(raw_data['orbit']['covariance']['data'])).astype(float)
-    # convert covariance matrix angle blocks from degrees to radians
-    cov_mat[3:6, :] *= np.pi/180
-    cov_mat[:, 3:6] *= np.pi/180
+    cov_exists = raw_data['orbit']['covariance'] is not None
+    cov_keys = raw_data['orbit']['covariance']['labels'] if cov_exists else []
+    if cov_exists:
+        cov_mat = (np.array(raw_data['orbit']['covariance']['data'])).astype(float)
+        # convert covariance matrix angle blocks from degrees to radians
+        cov_mat[3:6, :] *= np.pi/180
+        cov_mat[:, 3:6] *= np.pi/180
+    else:
+        print(f"WARNING: No covariance data available for {tdes}")
+        cov_mat = np.zeros((6, 6))
     # nongravitatinoal constants for target body
     nongrav_data = raw_data['orbit']['model_pars']
     hdr = [param['name'] for param in nongrav_data]