From 07e24556597637cde20b1cb495160a34329114fa Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Tue, 12 Apr 2022 12:26:16 -0400 Subject: [PATCH] start using commands-v2 --- agasc/supplement/magnitudes/mag_estimate.py | 15 +-- .../magnitudes/star_obs_catalogs.py | 124 +++++++++--------- agasc/tests/data/mag-stats.h5 | Bin 244368 -> 255328 bytes agasc/tests/test_obs_status.py | 13 -- 4 files changed, 67 insertions(+), 85 deletions(-) diff --git a/agasc/supplement/magnitudes/mag_estimate.py b/agasc/supplement/magnitudes/mag_estimate.py index 1925a9af..2d8f69fe 100644 --- a/agasc/supplement/magnitudes/mag_estimate.py +++ b/agasc/supplement/magnitudes/mag_estimate.py @@ -278,10 +278,9 @@ def get_telemetry(obs): :return: dict """ star_obs_catalogs.load() - dwell = star_obs_catalogs.DWELLS_NP[star_obs_catalogs.DWELLS_MAP[obs['mp_starcat_time']]] - star = get_star(obs['agasc_id'], date=dwell['tstart'], use_supplement=False) - start = dwell['tstart'] - stop = dwell['tstop'] + star = get_star(obs['agasc_id'], date=obs['obs_start'], use_supplement=False) + start = obs['obs_start'] + stop = obs['obs_stop'] slot = obs['slot'] logger.debug(f' Getting telemetry for AGASC ID={obs["agasc_id"]}, OBSID={obs["obsid"]}, ' f'mp_starcat_time={obs["mp_starcat_time"]}') @@ -638,9 +637,8 @@ def get_obs_stats(obs, telem=None): star_obs_catalogs.load() star = get_star(obs['agasc_id'], use_supplement=False) - dwell = star_obs_catalogs.DWELLS_NP[star_obs_catalogs.DWELLS_MAP[obs['mp_starcat_time']]] - start = dwell['tstart'] - stop = dwell['tstop'] + start = obs['obs_start'] + stop = obs['obs_stop'] stats = {k: obs[k] for k in ['agasc_id', 'obsid', 'slot', 'type', 'mp_starcat_time', 'timeline_id']} @@ -924,8 +922,7 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): star_obs_catalogs.load(tstop=tstop) # Get a table of every time the star has been observed - idx0, idx1 = star_obs_catalogs.STARS_OBS_MAP[agasc_id] - star_obs = star_obs_catalogs.STARS_OBS[idx0:idx1] + star_obs = star_obs_catalogs.STARS_OBS[star_obs_catalogs.STARS_OBS['agasc_id'] == agasc_id] if len(star_obs) > 1: star_obs = star_obs.loc['mp_starcat_time', sorted(star_obs['mp_starcat_time'])] diff --git a/agasc/supplement/magnitudes/star_obs_catalogs.py b/agasc/supplement/magnitudes/star_obs_catalogs.py index 5dcb724b..775b3db1 100644 --- a/agasc/supplement/magnitudes/star_obs_catalogs.py +++ b/agasc/supplement/magnitudes/star_obs_catalogs.py @@ -1,9 +1,8 @@ import os -import sys import numpy as np import tables -from astropy.table import Table +from astropy.table import Table, vstack, join from astropy import table from Ska.DBI import DBI @@ -11,20 +10,39 @@ from Chandra.Time import date2secs from chandra_aca.transform import yagzag_to_pixels -TIMELINES = None -STARCAT_CMDS = None -DWELLS_NP = None +_STARCAT_CMDS = None +_DWELLS = None """Array of dwells (mp_starcat_time, tstart, tstop, ra, dec, roll).""" -DWELLS_MAP = None -"""Time-index on DWELLS_NP. This is a dictionary mapping starcat times to an index in DWELLS_NP.""" - STARS_OBS = None """The table of star observations""" -STARS_OBS_INDICES = None -STARS_OBS_MAP = None -STARS_OBS_NP = None + +def get_star_observations(start=None, stop=None, obsid=None): + from kadi import commands + + with commands.conf.set_temp('commands_version', '2'): + catalogs = commands.get_starcats(start=start, stop=stop, obsid=obsid) + observations = Table(commands.get_observations(start=start, stop=stop, obsid=obsid)) + for cat in catalogs: + cat['obsid'] = cat.obsid + catalogs = vstack(catalogs, metadata_conflicts='silent') + catalogs = catalogs[np.in1d(catalogs['type'], ['BOT', 'GUI'])] + star_obs = join(observations, catalogs, keys=['obsid']) + star_obs.rename_columns(['id', 'starcat_date'], ['agasc_id', 'mp_starcat_time']) + star_obs['row'], star_obs['col'] = yagzag_to_pixels(star_obs['yang'], star_obs['zang']) + + # Add mag_aca_err column + + filename = os.path.join(os.environ['SKA'], 'data', 'agasc', 'proseco_agasc_1p7.h5') + with tables.open_file(filename) as h5: + agasc_ids = h5.root.data.col('AGASC_ID') + mag_errs = h5.root.data.col('MAG_ACA_ERR') * 0.01 + + tt = Table([agasc_ids, mag_errs], names=['agasc_id', 'mag_aca_err']) + star_obs = table.join(star_obs, tt, keys='agasc_id') + + return star_obs def _load_startcat_commands(tstop=None): @@ -49,15 +67,14 @@ def _load_startcat_commands(tstop=None): from kadi import events from kadi.commands import get_cmds - global STARCAT_CMDS - global DWELLS_NP - global DWELLS_MAP + global _STARCAT_CMDS + global _DWELLS - STARCAT_CMDS = Table(get_cmds('2003:001', type='MP_STARCAT'))['date', 'timeline_id'] - STARCAT_CMDS.convert_bytestring_to_unicode() - STARCAT_CMDS.rename_column('date', 'mp_starcat_time') + _STARCAT_CMDS = Table(get_cmds('2003:001', type='MP_STARCAT'))['date', 'timeline_id'] + _STARCAT_CMDS.convert_bytestring_to_unicode() + _STARCAT_CMDS.rename_column('date', 'mp_starcat_time') # print(len(starcat_cmds)) - STARCAT_CMDS = table.unique(STARCAT_CMDS, keys='mp_starcat_time') + _STARCAT_CMDS = table.unique(_STARCAT_CMDS, keys='mp_starcat_time') # print(len(starcat_cmds)) # Get all maneuvers and use this to create a table of NPM dwells @@ -73,34 +90,35 @@ def _load_startcat_commands(tstop=None): # Form dwells table. At this point it includes maneuvers that might not have # a star catalog with NPM dwell. That gets fixed in the next bit. dwells = Table(rows=rows, names=('obsid', 'manvr_start', 'start', 'stop', 'ra', 'dec', 'roll')) + dwells.convert_bytestring_to_unicode() # Clip starcat_cmds to be within time range of dwells - ok = (STARCAT_CMDS['mp_starcat_time'] < dwells['manvr_start'][-1]) - STARCAT_CMDS = STARCAT_CMDS[ok] + ok = (_STARCAT_CMDS['mp_starcat_time'] < dwells['manvr_start'][-1]) + _STARCAT_CMDS = _STARCAT_CMDS[ok] if tstop is not None: - STARCAT_CMDS = STARCAT_CMDS[STARCAT_CMDS['mp_starcat_time'] <= CxoTime(tstop).date] + _STARCAT_CMDS = _STARCAT_CMDS[_STARCAT_CMDS['mp_starcat_time'] <= CxoTime(tstop).date] # Now make a time-based association between every star catalog command and the # subsequent maneuver command. By the way that commands are built it is req'd # that a star catalog command be followed by a tlmsid=MANUVR command which # corresponds to `manvr_start` in the kadi manvrs events. Empirically this is # always less than 1000 sec delay (so allow up to 1200). - idxs = np.searchsorted(dwells['manvr_start'], STARCAT_CMDS['mp_starcat_time']) + idxs = np.searchsorted(dwells['manvr_start'], _STARCAT_CMDS['mp_starcat_time']) # Check that the association is correct. - dt = CxoTime(dwells['manvr_start'][idxs]) - CxoTime(STARCAT_CMDS['mp_starcat_time']) + dt = CxoTime(dwells['manvr_start'][idxs]) - CxoTime(_STARCAT_CMDS['mp_starcat_time']) ok = (dt.sec > 1) & (dt.sec < 1200) # print(f'WARNING: {np.sum(~ok)} entries where star-catalog/maneuver association is not right') # assert np.all((dt.sec > 1) & (dt.sec < 1200)) - STARCAT_CMDS = STARCAT_CMDS[ok] + _STARCAT_CMDS = _STARCAT_CMDS[ok] # Now re-select dwells to be only dwells with a star catalog and set the # corresponding command time for MP_STARCAT. This is the unique key to # allow going from a star catalog to the NPM dwell that uses the catalog. dwells = dwells[idxs][ok] - dwells['mp_starcat_time'] = STARCAT_CMDS['mp_starcat_time'] + dwells['mp_starcat_time'] = _STARCAT_CMDS['mp_starcat_time'] dwells.add_index('mp_starcat_time') # Kick out a few cases with None values for start or stop @@ -110,13 +128,7 @@ def _load_startcat_commands(tstop=None): # Now make a numpy array version for speed and a map (fast version of .loc) dwells['tstart'] = date2secs(dwells['start'].tolist()) dwells['tstop'] = date2secs(dwells['stop'].tolist()) - DWELLS_NP = dwells['mp_starcat_time', 'tstart', 'tstop', 'ra', 'dec', 'roll'].as_array() - if sys.version_info.minor == 6: - DWELLS_MAP = {DWELLS_NP['mp_starcat_time'][idx].decode('ascii'): idx - for idx in range(len(DWELLS_NP))} - else: - DWELLS_MAP = {DWELLS_NP['mp_starcat_time'][idx]: idx - for idx in range(len(DWELLS_NP))} + _DWELLS = dwells['mp_starcat_time', 'tstart', 'tstop', 'ra', 'dec', 'roll'] if 'DJANGO_SETTINGS_MODULE' in os.environ: del os.environ['DJANGO_SETTINGS_MODULE'] @@ -138,22 +150,19 @@ def _load_observed_catalogs(): if a star catalog entry *might* have been run. But not all approved load dirs were completed so we need to also cross-correlate with obsids from telemetry. """ - global TIMELINES global STARS_OBS - global STARS_OBS_INDICES global STARS_OBS_MAP - global STARS_OBS_NP # Timelines defines the mission planning `dir` of loads that were approved and # at least partially run. cmd_states_db3 = os.path.join(os.environ['SKA'], 'data', 'cmd_states', 'cmd_states.db3') with DBI(server=cmd_states_db3, dbi='sqlite') as timelines_db: - TIMELINES = Table(timelines_db.fetchall('select dir, id from timelines ' + timelines = Table(timelines_db.fetchall('select dir, id from timelines ' 'where datestart > "2003:007"')) - TIMELINES.rename_column('id', 'timeline_id') - TIMELINES.add_index('dir') - TIMELINES.add_index('timeline_id') + timelines.rename_column('id', 'timeline_id') + timelines.add_index('dir') + timelines.add_index('timeline_id') # Get table of every star catalog entry (a single star or fid light) that # was every planned. @@ -179,7 +188,7 @@ def _load_observed_catalogs(): ok = np.in1d(cat_entries['type'], ['GUI', 'BOT']) cat_entries = Table(cat_entries[ok]) cat_entries.rename_column('id', 'agasc_id') - cat_entries = table.join(cat_entries, TIMELINES, keys='dir') + cat_entries = table.join(cat_entries, timelines, keys='dir') # print(len(cat_entries)) # cat_entries[-5:].pprint(max_width=-1) @@ -187,7 +196,7 @@ def _load_observed_catalogs(): # joining `cat_entries` (which has a lot of candidates that might have been run) # with `starcat_cmds` that were actually run in a timeline. - ces = table.join(cat_entries, STARCAT_CMDS, keys=['mp_starcat_time', 'timeline_id']) + ces = table.join(cat_entries, _STARCAT_CMDS, keys=['mp_starcat_time', 'timeline_id']) # print(len(ces)) # print(ces[-20:]) # ces[ces['obsid'] == 19809] @@ -203,7 +212,7 @@ def _load_observed_catalogs(): ces = table.join(ces, tt, keys='agasc_id') # dropping catalog entries with not corresponding dwell # (these were usually dropped because start or stop was None) - ces = ces[np.in1d(ces['mp_starcat_time'], DWELLS_NP['mp_starcat_time'])] + ces = ces[np.in1d(ces['mp_starcat_time'], _DWELLS['mp_starcat_time'])] # print(ces[-5:]) # Make final table of observed guide star catalog entries `star_obs` @@ -216,28 +225,17 @@ def _load_observed_catalogs(): ces['col'] = col STARS_OBS = ces.group_by('agasc_id') + STARS_OBS = join(STARS_OBS, _DWELLS, keys=['mp_starcat_time']) + STARS_OBS.rename_columns(['tstart', 'tstop'], ['obs_start', 'obs_stop']) STARS_OBS.add_index('agasc_id') STARS_OBS.add_index('mp_starcat_time') - # Make a table that contains the indices into `star_obs` and gives the - # repeat count `n_obs` for each star. - - indices = STARS_OBS.groups.indices - rows = [] - for idx0, idx1 in zip(indices[:-1], indices[1:]): - agasc_id = STARS_OBS['agasc_id'][idx0] - rows.append((agasc_id, idx1 - idx0, idx0, idx1)) - STARS_OBS_INDICES = Table(rows=rows, names=['agasc_id', 'n_obs', 'idx0', 'idx1']) - STARS_OBS_INDICES.sort('n_obs') - STARS_OBS_INDICES.add_index('agasc_id') - # print(stars_obs_indices[:5]) - # print(stars_obs_indices[-10:]) - - STARS_OBS_MAP = {row['agasc_id']: (row['idx0'], row['idx1']) for row in STARS_OBS_INDICES} - STARS_OBS_NP = STARS_OBS.as_array() - -def load(tstop=None): - if STARCAT_CMDS is None: - _load_startcat_commands(tstop) - _load_observed_catalogs() +def load(tstop=None, commands_v2=False): + global STARS_OBS + if STARS_OBS is None: + if commands_v2: + STARS_OBS = get_star_observations(stop=tstop) + else: + _load_startcat_commands(tstop) + _load_observed_catalogs() diff --git a/agasc/tests/data/mag-stats.h5 b/agasc/tests/data/mag-stats.h5 index 8c8b98a7efb2adbd6d815907d75bb60de2a3e23c..63e25a3a665e6c0458f94cdc89b1d25836c827c0 100644 GIT binary patch delta 868 zcmYjOOH30{6n$@|(+q9ZezZ`d7K%tUL8q0P+Gv82+7|nxrN+%DMO>I*AW#xr5Ere{ zg=t!@8xt4SxX?rs$5{~E2um?h*tvuSM%}S+;m$krUSlUSbKW`UzBBi}Y5dv0y>Xxs z6ef4IQ(1|Qz5LJdo9UNA^}(AxV|0mTD7&8v|3jm)m4g*uu@k= zLA6nw$mr~~Wx3A_tHW+j-ePdChmg03vr`Lqrkh$neC>1S`y$^37b-UwD&_eKy2e#& zYJPa%eUKZob4(`W2{*dt%S?L8gd_Ahcb0o(`iieXJAe3OW5XDRkWDw+!QRDeZj@q| zG8Y>kMSUwUV8yXq5`0Kh;oY$yM3OZF1rmmDXT6|~&6x7toGHWBRfe=sgWp&_aE0L; zj(_uvKji-HEaQ8)->5I88LF8Y4Elqxo?_UFGh7;H_zU+L{4v@Z47AYQ!E^&JV!DAxG2OsL><&?+M1t|a8JM!aa82&A!&+*e6|75s K^H(?5e*Od1!q5@` delta 287 zcmaEGihsgaz6lyk6LxLX@?&MRm|V%~&dJCC0SXg0if?|vD!`(|0Tt0;hcI9&0-*c? z&dowxYK)U7h)PUvxXZkW(}4$~uy*oBCHcvLrec$4$cb(a<+;r>xzJL96RM4Yeey<8 z!_8fiEH;xnL=}|fq1M1?W`+s|kX2$#3=A?rIwiAc;zrf!7j`qVh>8IDC7HRYIhlE> z@tG-Logx!A%1=J0$TRt+wj5VqELe;og?X}|j@;x&igJ^CbrmMB5S3w6oBT0Nmh)2r fL_yK!jY_U;6AySa_r