""" Code to load the OSIRIS Level 1 data from the standard HDF4 / onyx files.
The methods open_level1_*() all perform the same basic function
of opening an OSIRIS Level 1 HDF4 data file, reading the data and
organizing it into an xarray dataset for convenience / usability.
There are optional parameters to load audit / extra data on some of the methods.
"""
from typing import Dict, List, NamedTuple
import sys
import os
from pathlib import Path
from functools import lru_cache
import xarray as xr
from numpy.testing import assert_equal, assert_array_equal
from osirisl1services.pyhdfadapter import HDF4VS
import numpy as np
import pandas as pd
from astropy.time import Time
from bisect import bisect
import warnings
import logging
# Use a large lru_cache if code is running in tests as it speeds things up. Otherwise, 4 is enough for most use cases.
if "pytest" in sys.modules:
logging.warning('Running in test mode. Using a limitless lru_cache for hdf4 files.')
LRU_CACHE_SIZE = None
else:
logging.info('Using an HDF4 cache size of 4.')
LRU_CACHE_SIZE = 4
MIN_RADIATION_HIT = 10
BLUE_RED_WAVE_PIXEL = 650
BLUE_RED_FWHM_PIXEL = 575
# If the attitude file orbit start time differs from the corresponding IR or spectrograph data start time by more than
# this time, load the previous orbit's data as well.
MAX_ATTITUDE_DATA_TIME_DELTA = 5.0 / 24.0 / 3600.0 # 5 seconds (mjd)
[docs]class Orbit(int):
""" L1 Orbit class with directory helpers.
:param int Orbit: the orbit number
"""
class L1FileParts(NamedTuple):
prefix: str
postfix: str
L1Dict = Dict[str, List[L1FileParts]]
l1dict: L1Dict = {'IR level1': [L1FileParts(prefix='si1_', postfix='.ir1')],
'IR level1b': [L1FileParts(prefix='si1b', postfix='.ir1')],
'attitude': [L1FileParts(prefix='so1_', postfix='.oat')],
'ecmwf': [L1FileParts(prefix='se1_', postfix='.ecm')],
'housekeeping': [
L1FileParts(prefix='sh0_', postfix='.hk1'),
L1FileParts(prefix='sh1_', postfix='.hk1')],
'radiation hit': [L1FileParts(prefix='sp1_', postfix='.txt')],
'spectrograph': [L1FileParts(prefix='ss1_', postfix='.os1')]}
def __new__(cls, value: int):
if value < 3713:
warnings.warn('Valid orbits start at 3713.')
return int.__new__(cls, int(value))
def as_hex(self):
return f'{self:04x}'
def l1_folder(self) -> str:
""" Return the base folder for the level 1 HDF4 / onyx files. """
cuts, folders = [9000, 19000], ['disk1', 'disk2', 'disk3']
folder = Path(os.environ['ODINORBITDIR'], folders[bisect(cuts, self)], f'{self // 100 * 100:05d}')
if not folder.exists():
folder = Path(os.environ['ODINORBITDIR'], f'{self // 100 * 100:05d}')
return folder
def _l1_filename(self, l1type: str) -> str:
""" Look up the filename for the given Level1 data type. """
prepost = Orbit.l1dict[l1type]
for (pre, post) in prepost:
filename = Path(self.l1_folder(), pre + self.as_hex() + post)
if filename.is_file():
return filename
else:
raise LookupError(f'{l1type} file not found in orbit {self} / {self:x}')
@property
def l1_attitude_filename(self):
return self._l1_filename('attitude')
@property
def l1_housekeeping_filename(self):
return self._l1_filename('housekeeping')
@property
def l1_spectrograph_filename(self):
return self._l1_filename('spectrograph')
@property
def l1_radiation_hit_filename(self):
return self._l1_filename('radiation hit')
@property
def l1_ir_filename(self):
return self._l1_filename('IR level1')
@property
def l1b_ir_filename(self):
return self._l1_filename('IR level1b')
@property
def l1_ecmwf_filename(self):
return self._l1_filename('ecmwf')
def _limit_to_scan(ds: xr.Dataset, scanno: int) -> xr.Dataset:
""" Reduce the returned dataset to a single scan's data """
att = open_level1_attitude(orbit=scanno // 1000)
start_mjd, end_mjd = (att.StartMJD.loc[scanno], att.EndMJD.loc[scanno])
reduced = ds.sel(mjd=slice(start_mjd, end_mjd))
if reduced.sizes['mjd'] == 0:
raise KeyError
return reduced
[docs]@lru_cache(maxsize=LRU_CACHE_SIZE)
def open_level1_attitude(orbit: int, scanno: int = None, load_pointing=False, load_ssc_los=False) -> xr.Dataset:
""" Load the OSIRIS Level 1 attitude data for an orbit (or scan). Optinally load the pointing and SSC's LOS data.
:param Orbit orbit: The orbit class handling the data retrieval, specified by orbit number
:param int scanno: The scan number to limit the data to, default is None
:param bool load_pointing: If True, pointing data is included in the dataset
:param bool load_ssc_los: If True, SSC's line of sight data is included in the dataset
The returned data structure contains:
<xarray.Dataset>
Dimensions: (Q: 4, ScanNumber: 66, mjd: 5464, xyz: 3)
Coordinates:
* mjd (mjd) float64 5.239e+04 5.239e+04 5.239e+04 ...
* xyz (xyz) <U1 'x' 'y' 'z'
* ScanNumber (ScanNumber) int64 6432000 6432001 6432002 6432003 ...
Data variables:
mjdorbitstate (mjd) float64 5.239e+04 5.239e+04 5.239e+04 ...
** the mjd used to compute the orbit state parameters (position, velocity) **
stw (mjd) int64 598896805 598896822 598896839 ...
** satellite time word **
position (mjd, xyz) float64 -4.252e+06 5.54e+06 1.196e+03 ...
** satellite position in ECI coordinate system **
velocity (mjd, xyz) float64 823.0 615.0 7.485e+03 828.0 ...
** satellite velocity in ECI coordinate system **
version (mjd) float64 17.1 17.1 17.1 17.1 17.1 17.1 17.1 ...
** The current version of satellite software? **
StartMJD (ScanNumber) float64 5.239e+04 5.239e+04 5.239e+04 ...
** Start time of each scan in modified Julian date **
EndMJD (ScanNumber) float64 5.239e+04 5.239e+04 5.239e+04 ...
** End time of each scan in modified Julian date **
ACSState (ScanNumber) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
** The current state of the ACS system. 1. limb scanning 2. limb staring 3. other **
MinAltitude (ScanNumber) float64 14.43 14.41 14.41 14.4 14.4 ...
MaxAltitude (ScanNumber) float64 76.24 76.24 76.26 76.26 76.31 ...
** min / max altitude of the control frame, not either of the instruments
ScanRate (ScanNumber) float64 0.7095 0.7098 0.7099 0.71 ...
** km / sec
Direction (ScanNumber) int64 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 ...
** 1. Up 2. Down **
If the pointing data is loaded, the returned data structure also contains:
Coordinates:
* 1ijk (1ijk) <U1 '1' 'i' 'j' 'k'
Data variables:
Qtarget (mjd, 1ijk) float64 0.7487 0.335 -0.5648 0.0901 0.7488 ...
Qachieved (mjd, 1ijk) float64 -0.7487 -0.335 0.5649 -0.0902 -0.7488 ...
Qerror (mjd, 1ijk) float64 15.0 5.0 3.0 0.0 15.0 6.0 3.0 0.0 ...
If the SSC line of sight data is loaded, the returned data structure also contains:
Data variables:
osiris_los_latitude (mjd) float64 23.72 23.78 23.82 23.86 23.91 23.95 ...
** osiris latitude data computed by SSC (differs from level1services) **
osiris_los_longitude (mjd) float64 134.9 134.9 134.9 134.9 134.9 134.9 ...
** osiris longitude data computed by SSC (differs from level1 services) **
osiris_los_height (mjd) float64 7.989 8.956 9.889 10.76 11.59 12.4 ...
** osiris altitude data computed by SSC (differs from level1 services) **
For example:
>>> attitude = open_level1_attitude(orbit=6432)
"""
with HDF4VS(Orbit(orbit).l1_attitude_filename) as vs:
sa = vs['Spacecraft Attitude']
st = vs['Scan Times']
ds = xr.Dataset({'mjdorbitstate': (['mjd'], sa['mjdorbitstate']),
'stw': (['mjd'], sa['stw']),
'position': (['mjd', 'xyz'], sa['position']),
'velocity': (['mjd', 'xyz'], sa['velocity']),
'version': (['mjd'], sa['version']),
'StartMJD': (['ScanNumber'], st['StartMJD']),
'EndMJD': (['ScanNumber'], st['EndMJD']),
'ACSState': (['ScanNumber'], st['ACSState']),
# 'ACSConfig': (['ScanNumber'], st['ACSConfig']), # removed (not properly implemented)
'MinAltitude': (['ScanNumber'], st['MinAltitude']),
'MaxAltitude': (['ScanNumber'], st['MaxAltitude']),
'ScanRate': (['ScanNumber'], st['ScanRate']),
'Direction': (['ScanNumber'], st['Direction'])},
coords={'mjd': sa['mjd'],
'xyz': ['x', 'y', 'z'],
'ScanNumber': st['ScanNumber']})
if load_pointing:
pointing = xr.Dataset({'Qtarget': (['mjd', '1ijk'], sa['Qtarget']),
'Qachieved': (['mjd', '1ijk'], sa['Qachieved']),
'Qerror': (['mjd', '1ijk'], sa['Qerror']),
},
coords={'mjd': sa['mjd'],
'1ijk': ['1', 'i', 'j', 'k']})
ds = ds.merge(pointing)
if load_ssc_los:
line_of_sight = xr.Dataset({'osiris_los_latitude': (['mjd'], sa['osiris_los_latitude']),
'osiris_los_longitude': (['mjd'], sa['osiris_los_longitude']),
'osiris_los_height': (['mjd'], sa['osiris_los_height']),
# smr data not used by our group
# 'smr_los_latitude': (['mjd'], sa['smr_los_latitude']),
# 'smr_los_longitude': (['mjd'], sa['smr_los_longitude']),
# 'smr_los_height': (['mjd'], sa['smr_los_height']),
},
coords={'mjd': sa['mjd']})
ds = ds.merge(line_of_sight)
# Reduce the amount of data returned to a single scan, if given
# As the attitude file is already open, this method is slightly different than _limit_to_scan(ds, orbit, scanno)
if scanno is not None:
try:
start_mjd, end_mjd = (ds.StartMJD.loc[scanno], ds.EndMJD.loc[scanno])
except KeyError as e:
print(f'Scan {scanno} is not in Orbit {orbit}')
raise e
# return the attitude information for the scan with a small buffer on either side to compensate for the fact
# that the attitude and data mjds to not align perfectly
return ds.sel(mjd=slice(start_mjd - MAX_ATTITUDE_DATA_TIME_DELTA,
end_mjd + MAX_ATTITUDE_DATA_TIME_DELTA), ScanNumber=scanno)
return ds
[docs]@lru_cache(maxsize=LRU_CACHE_SIZE)
def open_level1_housekeeping(orbit: int) -> xr.Dataset:
""" Load the OSIRIS Level 1 housekeeping data for an orbit.
:param Orbit orbit: The orbit id to load data from
:param int scanno: The scan number to limit data to, if given
The returned data structure contains:
<xarray.Dataset>
Dimensions: (dettempmjd: 1940, mjd: 90, statmjd: 39)
Coordinates:
* mjd (mjd) float64 5.239e+04 5.239e+04 5.239e+04 ...
* statmjd (statmjd) float64 5.239e+04 5.239e+04 5.239e+04 ...
* dettempmjd (dettempmjd) float64 5.239e+04 5.239e+04 ...
Data variables:
stw (mjd) int64 598897088 598898128 598899168 ...
plus_28v (mjd) float64 30.14 30.16 30.11 30.16 30.16 ...
milliamps (mjd) float64 569.3 565.0 575.7 568.2 501.1 ...
plus_5v (mjd) float64 4.947 4.928 4.91 4.91 4.928 4.91 ...
plus_12v (mjd) float64 11.81 11.85 11.85 11.81 11.85 ...
plus_24v (mjd) float64 23.53 23.72 23.53 23.53 23.53 ...
plus_36v (mjd) float64 35.6 35.47 35.6 35.33 35.47 35.6 ...
neg_5v (mjd) float64 -5.002 -4.984 -4.984 -5.002 -5.002 ...
neg_12v (mjd) float64 -12.08 -12.08 -11.99 -12.03 -11.99 ...
gnd_mux1 (mjd) float64 0.7912 0.8029 0.7971 0.7853 0.7971 ...
gnd_mux2 (mjd) float64 0.0119 0.01776 0.0119 0.0119 ...
CPU (mjd) float64 49.53 48.79 49.16 49.16 48.79 ...
PS1 (mjd) float64 65.17 65.83 67.2 65.17 65.83 65.83 ...
PS2 (mjd) float64 67.92 68.65 67.92 67.2 67.2 68.65 ...
IRStrap (mjd) float64 -13.5 -13.34 -13.5 -13.5 -13.5 ...
OSStrap (mjd) float64 -18.43 -18.36 -18.43 -18.36 -18.36 ...
IFPlate (mjd) float64 26.85 26.85 26.85 26.68 27.03 ...
Optics (mjd) float64 22.83 22.83 23.15 22.83 22.99 ...
IR1 (mjd) float64 -13.17 -13.09 -12.84 -13.0 -13.09 ...
IR2 (mjd) float64 -17.29 -17.51 -17.35 -17.35 -17.43 ...
IR3 (mjd) float64 -16.4 -16.4 -16.31 -16.24 -16.4 ...
CCD (mjd) float64 -12.54 -12.54 -12.62 -12.46 -12.46 ...
RadiatorPlate (mjd) float64 -30.76 -30.69 -30.69 -30.69 -30.69 ...
OS_TecIsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
OS_TecDutyCycle (mjd) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
IR_TecIsOn (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
IR_TecDutyCycle (mjd) float64 74.74 74.74 74.74 74.74 74.74 ...
CCD_HtrIsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
IR1_HtrIsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
IR2_HtrIsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
IR3_HtrIsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
CCD_HtrSetPoint (mjd) int64 14179 14179 14179 14179 14179 14179 ...
IR1_HtrSetPoint (mjd) int64 15360 15360 15360 15360 15360 15360 ...
IR2_HtrSetPoint (mjd) int64 15360 15360 15360 15360 15360 15360 ...
IR3_HtrSetPoint (mjd) int64 15360 15360 15360 15360 15360 15360 ...
IR_State (mjd) int64 10 9 10 10 10 10 10 10 10 10 10 10 ...
OS_State (mjd) int64 10 10 10 10 10 10 10 10 10 10 10 10 ...
IR_SciProg (mjd) int64 34062 34062 34062 34062 34062 34062 ...
OS_SciProg (mjd) int64 21 21 21 21 21 21 21 21 21 21 21 21 ...
SPM_IsBusy (mjd) int64 0 0 1 1 0 1 1 0 0 0 1 0 0 0 1 1 1 0 ...
TargetIndex (mjd) int64 57879 58935 59727 57879 59463 57615 ...
LastError (mjd) int64 286 286 286 286 286 286 286 286 286 ...
NumErrors (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
OSShutterPos (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
IRShutterPos (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
IRShutterDriverIsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
OSShutterDriverIsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
IRLamp1IsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
IRLamp2IsOn (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
SEUTriggerOccured (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
SEUDetectorIsActive (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
SEUCounter (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
FifoStatus (mjd) int64 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 ...
FifoOverrun (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
FormatsAreLost (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
MassMemoryIsAvailable (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
MassMemoryCounter (mjd) int64 6903 19004 18865 18648 18467 18276 ...
SciWordsPerFormat (mjd) int64 156 156 156 156 156 156 156 156 156 ...
DataIsPendingInSciBuf (mjd) int64 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
UpperFlashOk (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
LowerFlashOk (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
FlashLoadOk (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
InOrbitEnabled (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
VerboseMessagesEnabled (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
TraceEnabled (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
UsingACDC1 (mjd) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
UsingACDC2 (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
BootSource (mjd) int64 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 ...
fswrev (mjd) int64 50 50 50 50 50 50 50 50 50 50 50 50 ...
ROM_Checksum (mjd) int64 4164 16452 16452 16452 16452 16452 ...
LowerFlash_Checksum (mjd) int64 227357 247837 202781 100381 145437 ...
UpperFlash_Checksum (mjd) int64 177071 164783 156591 144303 132015 ...
osiris_stw (mjd) int64 1209092 1210132 1211172 1212212 ...
LastCommand (mjd) int64 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 ...
NumCommandWordsPending (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
NumBadCommands (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
NumMarginalCommands (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
NumLostCommandwords (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
NumTotalCommands (mjd) int64 49 49 49 49 49 49 49 49 49 49 49 49 ...
statstw (statmjd) int64 598897404 598899804 598902204 ...
statsnumpoints (statmjd) int64 2400 2400 2400 2400 2400 2400 ...
plus5v_minval (statmjd) float64 4.883 4.852 4.852 4.881 4.853 ...
plus5v_maxval (statmjd) float64 5.041 5.005 5.004 5.022 5.022 ...
plus5v_average (statmjd) float64 4.934 4.933 4.933 4.934 4.934 ...
plus12v_minval (statmjd) float64 11.73 11.73 11.72 11.73 11.72 ...
plus12v_maxval (statmjd) float64 12.04 12.01 12.0 11.99 12.01 ...
plus12v_average (statmjd) float64 11.85 11.85 11.85 11.85 11.85 ...
plus24v_minval (statmjd) float64 23.34 23.34 23.36 23.43 23.25 ...
plus24v_maxval (statmjd) float64 24.08 24.1 24.1 24.01 24.19 ...
plus24v_average (statmjd) float64 23.66 23.65 23.66 23.66 23.66 ...
plus36v_minval (statmjd) float64 34.95 34.95 34.95 35.09 34.92 ...
plus36v_maxval (statmjd) float64 36.14 36.04 36.01 36.14 36.01 ...
plus36v_average (statmjd) float64 35.5 35.5 35.5 35.5 35.5 35.51 ...
neg5v_minval (statmjd) float64 -5.041 -5.037 -5.055 -5.055 ...
neg5v_maxval (statmjd) float64 -4.913 -4.927 -4.908 -4.89 ...
neg5v_average (statmjd) float64 -4.979 -4.979 -4.979 -4.979 ...
neg12v_minval (statmjd) float64 -12.2 -12.14 -12.14 -12.17 ...
neg12v_maxval (statmjd) float64 -11.82 -11.87 -11.84 -11.87 ...
neg12v_average (statmjd) float64 -12.01 -12.01 -12.01 -12.01 ...
gnd_mux1_minval (statmjd) float64 0.7703 0.7703 0.7659 0.7736 ...
gnd_mux1_maxval (statmjd) float64 0.822 0.8322 0.8322 0.8234 ...
gnd_mux1_average (statmjd) float64 0.7982 0.7985 0.7993 0.7996 ...
gnd_mux2_minval (statmjd) float64 -0.01044 -0.001648 -0.01044 ...
gnd_mux2_maxval (statmjd) float64 0.03534 0.03021 0.0412 0.03534 ...
gnd_mux2_average (statmjd) float64 0.0141 0.0141 0.0141 0.0141 ...
CPU_minval (statmjd) float64 47.6 47.86 47.77 48.03 48.06 ...
CPU_maxval (statmjd) float64 50.29 50.29 51.07 50.68 50.36 ...
CPU_average (statmjd) float64 49.09 49.16 49.2 49.27 49.3 ...
PS1_minval (statmjd) float64 63.23 63.08 63.38 63.54 63.38 ...
PS1_maxval (statmjd) float64 68.1 68.01 69.21 68.65 67.92 ...
PS1_average (statmjd) float64 65.83 65.83 65.83 65.83 65.83 ...
PS2_minval (statmjd) float64 66.21 65.79 65.29 65.79 65.12 ...
PS2_maxval (statmjd) float64 73.05 71.37 71.42 70.97 71.27 ...
PS2_average (statmjd) float64 68.65 68.7 68.65 68.7 68.7 ...
IRStrap_minval (statmjd) float64 -13.77 -13.73 -13.83 -13.77 ...
IRStrap_maxval (statmjd) float64 -13.01 -13.09 -13.17 -13.17 ...
IRStrap_average (statmjd) float64 -13.42 -13.42 -13.43 -13.43 ...
OSStrap_minval (statmjd) float64 -18.63 -18.61 -18.64 -18.68 ...
OSStrap_maxval (statmjd) float64 -18.04 -18.1 -18.13 -18.13 ...
OSStrap_average (statmjd) float64 -18.36 -18.37 -18.37 -18.36 ...
IFPlate_minval (statmjd) float64 26.23 26.32 26.19 26.49 26.23 ...
IFPlate_maxval (statmjd) float64 27.65 27.67 27.44 27.57 27.57 ...
IFPlate_average (statmjd) float64 26.93 26.93 26.93 26.94 26.95 ...
Optics_minval (statmjd) float64 22.38 22.46 22.34 22.34 22.5 ...
Optics_maxval (statmjd) float64 23.47 23.49 23.5 23.63 23.47 ...
Optics_average (statmjd) float64 22.95 22.94 22.93 22.93 22.93 ...
IR1_minval (statmjd) float64 -13.35 -13.35 -13.42 -13.43 ...
IR1_maxval (statmjd) float64 -12.74 -12.72 -12.74 -12.75 ...
IR1_average (statmjd) float64 -13.02 -13.05 -13.07 -13.08 ...
IR2_minval (statmjd) float64 -17.69 -17.6 -17.72 -17.67 ...
IR2_maxval (statmjd) float64 -17.05 -17.1 -17.12 -17.12 ...
IR2_average (statmjd) float64 -17.35 -17.37 -17.39 -17.41 ...
IR3_minval (statmjd) float64 -16.57 -16.63 -16.65 -16.63 ...
IR3_maxval (statmjd) float64 -16.03 -15.99 -15.99 -16.07 ...
IR3_average (statmjd) float64 -16.29 -16.31 -16.33 -16.34 ...
CCD_minval (statmjd) float64 -12.79 -12.81 -12.89 -12.89 ...
CCD_maxval (statmjd) float64 -12.29 -12.22 -12.29 -12.29 ...
CCD_average (statmjd) float64 -12.53 -12.55 -12.57 -12.57 ...
dettempstw (dettempmjd) int64 598896850 598896897 598896945 ...
dettempCCD_T (dettempmjd) float64 -12.52 -12.52 -12.53 -12.53 ...
dettempIR1_T (dettempmjd) float64 -13.01 -13.02 -13.02 -13.01 ...
dettempIR2_T (dettempmjd) float64 -17.33 -17.34 -17.35 -17.34 ...
dettempIR3_T (dettempmjd) float64 -16.26 -16.27 -16.27 -16.3 ...
For example:
>>> housekeeping = open_level1_housekeeping(orbit=6432)
"""
with HDF4VS(Orbit(orbit).l1_housekeeping_filename) as vs:
v = vs['Voltages']
t = vs['Temperatures']
f = vs['HK Flags']
vstats = vs['Voltage Stats']
tstats = vs['Temperature Stats']
dt = vs['Det.Temp (Run. Avg)']
v = xr.Dataset({'stw': (['mjd'], v['stw']),
'plus_28v': (['mjd'], v['plus_28v']),
'milliamps': (['mjd'], v['milliamps']),
'plus_5v': (['mjd'], v['plus_5v']),
'plus_12v': (['mjd'], v['plus_12v']),
'plus_24v': (['mjd'], v['plus_24v']),
'plus_36v': (['mjd'], v['plus_36v']),
'neg_5v': (['mjd'], v['neg_5v']),
'neg_12v': (['mjd'], v['neg_12v']),
'gnd_mux1': (['mjd'], v['gnd_mux1']),
'gnd_mux2': (['mjd'], v['gnd_mux2'])},
coords={'mjd': v['mjd']})
t = xr.Dataset({'stw': (['mjd'], t['stw']),
'CPU': (['mjd'], t['CPU']),
'PS1': (['mjd'], t['PS1']),
'PS2': (['mjd'], t['PS2']),
'IRStrap': (['mjd'], t['IRStrap']),
'OSStrap': (['mjd'], t['OSStrap']),
'IFPlate': (['mjd'], t['IFPlate']),
'Optics': (['mjd'], t['Optics']),
'IR1': (['mjd'], t['IR1']),
'IR2': (['mjd'], t['IR2']),
'IR3': (['mjd'], t['IR3']),
'CCD': (['mjd'], t['CCD']),
'RadiatorPlate': (['mjd'], t['RadiatorPlate']),
'OS_TecIsOn': (['mjd'], t['OS_TecIsOn']),
'OS_TecDutyCycle': (['mjd'], t['OS_TecDutyCycle']),
'IR_TecIsOn': (['mjd'], t['IR_TecIsOn']),
'IR_TecDutyCycle': (['mjd'], t['IR_TecDutyCycle']),
'CCD_HtrIsOn': (['mjd'], t['CCD_HtrIsOn']),
'IR1_HtrIsOn': (['mjd'], t['IR1_HtrIsOn']),
'IR2_HtrIsOn': (['mjd'], t['IR2_HtrIsOn']),
'IR3_HtrIsOn': (['mjd'], t['IR3_HtrIsOn']),
'CCD_HtrSetPoint': (['mjd'], t['CCD_HtrSetPoint']),
'IR1_HtrSetPoint': (['mjd'], t['IR1_HtrSetPoint']),
'IR2_HtrSetPoint': (['mjd'], t['IR2_HtrSetPoint']),
'IR3_HtrSetPoint': (['mjd'], t['IR3_HtrSetPoint'])},
coords={'mjd': t['mjd']})
f = xr.Dataset({'stw': (['mjd'], f['stw']),
'IR_State': (['mjd'], f['IR_State']),
'OS_State': (['mjd'], f['OS_State']),
'IR_SciProg': (['mjd'], f['IR_SciProg']),
'OS_SciProg': (['mjd'], f['OS_SciProg']),
'SPM_IsBusy': (['mjd'], f['SPM_IsBusy']),
'TargetIndex': (['mjd'], f['TargetIndex']),
'LastError': (['mjd'], f['LastError']),
'NumErrors': (['mjd'], f['NumErrors']),
'OSShutterPos': (['mjd'], f['OSShutterPos']),
'IRShutterPos': (['mjd'], f['IRShutterPos']),
'IRShutterDriverIsOn': (['mjd'], f['IRShutterDriverIsOn']),
'OSShutterDriverIsOn': (['mjd'], f['OSShutterDriverIsOn']),
'IRLamp1IsOn': (['mjd'], f['IRLamp1IsOn']),
'IRLamp2IsOn': (['mjd'], f['IRLamp2IsOn']),
'SEUTriggerOccured': (['mjd'], f['SEUTriggerOccured']),
'SEUDetectorIsActive': (['mjd'], f['SEUDetectorIsActive']),
'SEUCounter': (['mjd'], f['SEUCounter']),
'FifoStatus': (['mjd'], f['FifoStatus']),
'FifoOverrun': (['mjd'], f['FifoOverrun']),
'FormatsAreLost': (['mjd'], f['FormatsAreLost']),
'MassMemoryIsAvailable': (['mjd'], f['MassMemoryIsAvailable']),
'MassMemoryCounter': (['mjd'], f['MassMemoryCounter']),
'SciWordsPerFormat': (['mjd'], f['SciWordsPerFormat']),
'DataIsPendingInSciBuf': (['mjd'], f['DataIsPendingInSciBuf']),
'UpperFlashOk': (['mjd'], f['UpperFlashOk']),
'LowerFlashOk': (['mjd'], f['LowerFlashOk']),
'FlashLoadOk': (['mjd'], f['FlashLoadOk']),
'InOrbitEnabled': (['mjd'], f['InOrbitEnabled']),
'VerboseMessagesEnabled': (['mjd'], f['VerboseMessagesEnabled']),
'TraceEnabled': (['mjd'], f['TraceEnabled']),
'UsingACDC1': (['mjd'], f['UsingACDC1']),
'UsingACDC2': (['mjd'], f['UsingACDC2']),
'BootSource': (['mjd'], f['BootSource']),
'fswrev': (['mjd'], f['fswrev']),
'ROM_Checksum': (['mjd'], f['ROM_Checksum']),
'LowerFlash_Checksum': (['mjd'], f['LowerFlash_Checksum']),
'UpperFlash_Checksum': (['mjd'], f['UpperFlash_Checksum']),
'osiris_stw': (['mjd'], f['osiris_stw']),
'LastCommand': (['mjd'], f['LastCommand']),
'NumCommandWordsPending': (['mjd'], f['NumCommandWordsPending']),
'NumBadCommands': (['mjd'], f['NumBadCommands']),
'NumMarginalCommands': (['mjd'], f['NumMarginalCommands']),
'NumLostCommandwords': (['mjd'], f['NumLostCommandwords']),
'NumTotalCommands': (['mjd'], f['NumTotalCommands'])},
coords={'mjd': f['mjd']})
vstats = xr.Dataset({'statstw': (['statmjd'], vstats['stw']),
'statsnumpoints': (['statmjd'], vstats['numpoints']),
'plus5v_minval': (['statmjd'], vstats['plus5v_minval']),
'plus5v_maxval': (['statmjd'], vstats['plus5v_maxval']),
'plus5v_average': (['statmjd'], vstats['plus5v_average']),
'plus12v_minval': (['statmjd'], vstats['plus12v_minval']),
'plus12v_maxval': (['statmjd'], vstats['plus12v_maxval']),
'plus12v_average': (['statmjd'], vstats['plus12v_average']),
'plus24v_minval': (['statmjd'], vstats['plus24v_minval']),
'plus24v_maxval': (['statmjd'], vstats['plus24v_maxval']),
'plus24v_average': (['statmjd'], vstats['plus24v_average']),
'plus36v_minval': (['statmjd'], vstats['plus36v_minval']),
'plus36v_maxval': (['statmjd'], vstats['plus36v_maxval']),
'plus36v_average': (['statmjd'], vstats['plus36v_average']),
'neg5v_minval': (['statmjd'], vstats['neg5v_minval']),
'neg5v_maxval': (['statmjd'], vstats['neg5v_maxval']),
'neg5v_average': (['statmjd'], vstats['neg5v_average']),
'neg12v_minval': (['statmjd'], vstats['neg12v_minval']),
'neg12v_maxval': (['statmjd'], vstats['neg12v_maxval']),
'neg12v_average': (['statmjd'], vstats['neg12v_average']),
'gnd_mux1_minval': (['statmjd'], vstats['gnd_mux1_minval']),
'gnd_mux1_maxval': (['statmjd'], vstats['gnd_mux1_maxval']),
'gnd_mux1_average': (['statmjd'], vstats['gnd_mux1_average']),
'gnd_mux2_minval': (['statmjd'], vstats['gnd_mux2_minval']),
'gnd_mux2_maxval': (['statmjd'], vstats['gnd_mux2_maxval']),
'gnd_mux2_average': (['statmjd'], vstats['gnd_mux2_average'])},
coords={'statmjd': vstats['mjd']})
tstats = xr.Dataset({'statstw': (['statmjd'], tstats['stw']),
'statsnumpoints': (['statmjd'], tstats['numpoints']),
'CPU_minval': (['statmjd'], tstats['CPU_minval']),
'CPU_maxval': (['statmjd'], tstats['CPU_maxval']),
'CPU_average': (['statmjd'], tstats['CPU_average']),
'PS1_minval': (['statmjd'], tstats['PS1_minval']),
'PS1_maxval': (['statmjd'], tstats['PS1_maxval']),
'PS1_average': (['statmjd'], tstats['PS1_average']),
'PS2_minval': (['statmjd'], tstats['PS2_minval']),
'PS2_maxval': (['statmjd'], tstats['PS2_maxval']),
'PS2_average': (['statmjd'], tstats['PS2_average']),
'IRStrap_minval': (['statmjd'], tstats['IRStrap_minval']),
'IRStrap_maxval': (['statmjd'], tstats['IRStrap_maxval']),
'IRStrap_average': (['statmjd'], tstats['IRStrap_average']),
'OSStrap_minval': (['statmjd'], tstats['OSStrap_minval']),
'OSStrap_maxval': (['statmjd'], tstats['OSStrap_maxval']),
'OSStrap_average': (['statmjd'], tstats['OSStrap_average']),
'IFPlate_minval': (['statmjd'], tstats['IFPlate_minval']),
'IFPlate_maxval': (['statmjd'], tstats['IFPlate_maxval']),
'IFPlate_average': (['statmjd'], tstats['IFPlate_average']),
'Optics_minval': (['statmjd'], tstats['Optics_minval']),
'Optics_maxval': (['statmjd'], tstats['Optics_maxval']),
'Optics_average': (['statmjd'], tstats['Optics_average']),
'IR1_minval': (['statmjd'], tstats['IR1_minval']),
'IR1_maxval': (['statmjd'], tstats['IR1_maxval']),
'IR1_average': (['statmjd'], tstats['IR1_average']),
'IR2_minval': (['statmjd'], tstats['IR2_minval']),
'IR2_maxval': (['statmjd'], tstats['IR2_maxval']),
'IR2_average': (['statmjd'], tstats['IR2_average']),
'IR3_minval': (['statmjd'], tstats['IR3_minval']),
'IR3_maxval': (['statmjd'], tstats['IR3_maxval']),
'IR3_average': (['statmjd'], tstats['IR3_average']),
'CCD_minval': (['statmjd'], tstats['CCD_minval']),
'CCD_maxval': (['statmjd'], tstats['CCD_maxval']),
'CCD_average': (['statmjd'], tstats['CCD_average'])},
coords={'statmjd': tstats['mjd']})
dt = xr.Dataset({'dettempstw': (['dettempmjd'], dt['stw']),
'dettempCCD_T': (['dettempmjd'], dt['CCD_T']),
'dettempIR1_T': (['dettempmjd'], dt['IR1_T']),
'dettempIR2_T': (['dettempmjd'], dt['IR2_T']),
'dettempIR3_T': (['dettempmjd'], dt['IR3_T'])},
coords={'dettempmjd': dt['mjd']})
ds = xr.merge([v, t, f, vstats, tstats, dt])
return ds
@lru_cache(maxsize=LRU_CACHE_SIZE)
def _open_level1_spectrograph_hdf4(o: int, load_audit) -> xr.Dataset:
with HDF4VS(Orbit(o).l1_spectrograph_filename) as vs:
l1 = vs['OS Level 1']
data = [vs['R8:R4:01353']['Array'][i[0]] for i in l1['data']]
assert_equal([vs['R8:R4:01353']['crosscheck'][i[0]] for i in l1['data']], [d[3] for d in l1['data']])
error = [vs['R8:R4:01353']['Array'][i[0]] for i in l1['error']]
assert_equal([vs['R8:R4:01353']['crosscheck'][i[0]] for i in l1['error']], [d[3] for d in l1['error']])
flags = [vs['UI1:UI1:01353']['Array'][i[0]] for i in l1['flags']]
assert_equal([vs['UI1:UI1:01353']['crosscheck'][i[0]] for i in l1['flags']], [d[3] for d in l1['flags']])
ds = xr.Dataset({'stw': (['mjd'], l1['stw']),
'exposureTime': (['mjd'], l1['exposureTime']),
'temperature': (['mjd'], l1['temperature']),
'tempavg': (['mjd'], l1['tempavg']),
'opticstemp': (['mjd'], l1['opticstemp']),
'straptemp': (['mjd'], l1['straptemp']),
'mode': (['mjd'], l1['mode']),
'scienceprog': (['mjd'], l1['scienceprog']),
'roe': (['mjd'], l1['roe']),
'shuttermode': (['mjd'], l1['shuttermode']),
'spm_baserow': (['mjd'], l1['spm_baserow']),
'spm_numrows': (['mjd'], l1['spm_numrows']),
'spm_processingMode': (['mjd'], l1['spm_processingMode']),
'targetIndex': (['mjd'], l1['targetIndex']),
'exceptions': (['mjd'], l1['exceptions']),
'processingflags': (['mjd'], l1['processingflags']),
'numcolumns': (['mjd'], l1['numcolumns']),
'numrows': (['mjd'], l1['numrows']),
'data': (['mjd', 'pixel'], data),
'error': (['mjd', 'pixel'], error),
'flags': (['mjd', 'pixel'], flags)},
coords={'mjd': l1['mjd'],
'pixel': range(1353)})
if load_audit:
audit = xr.Dataset({'audit': (['mjd', 'store'], l1['audit'])},
coords={'mjd': l1['mjd'],
'store': [0, 1, 2, 3]})
ds = ds.merge(audit)
return ds
@lru_cache()
def _nominal_spectrograph_wavelengths_and_fwhm() -> xr.Dataset:
blue_wave = np.polyval([-1.3557105e-8, 1.8603574e-5, 0.38731255, 274.96731], range(BLUE_RED_WAVE_PIXEL))
red_wave = np.polyval([4.0389e-6, 0.3949776, 273.3871444], range(BLUE_RED_WAVE_PIXEL, 1353))
wavelength = xr.DataArray(np.concatenate([blue_wave, red_wave]), dims='pixel')
blue_fwhm = np.polyval([1.64967016e-8, -1.54021711e-5, 3.60834321e-3, 7.76274366e-1], range(BLUE_RED_FWHM_PIXEL))
red_fwhm = np.polyval([3.031138e-9, -1.081015e-5, 1.173076e-2, -2.76315896], range(BLUE_RED_FWHM_PIXEL, 1353))
fwhm = xr.DataArray(np.concatenate([blue_fwhm, red_fwhm]), dims='pixel')
return wavelength, fwhm
@lru_cache()
def _load_radiation_hit_data(orbit: int) -> xr.Dataset:
radiation_hits = pd.read_table(
Orbit(orbit).l1_radiation_hit_filename,
header=None, usecols=[0, 1, 2, 3], delim_whitespace=True,
names=['scannumber', 'stw', 'pixel', 'magnitude'],
dtype={'scannumber': int, 'stw': int, 'pixel': int, 'magnitude': float})
return radiation_hits
@lru_cache()
def _load_angular_rasa_data() -> xr.Dataset:
rsas_file = Path(os.environ['ODINFLIGHTDIR'], 'rsas_correction.nc')
rsas_data = xr.open_dataset(rsas_file)
rsas_data['mjd'] = ('time', Time(rsas_data.time).mjd + 0.5) # Shift to middle of day since rsas values are binned
return -rsas_data.swap_dims({'time': 'mjd'}).rsas_shift_rad_total.dropna(dim='mjd')
[docs]@lru_cache(maxsize=LRU_CACHE_SIZE)
def open_level1_spectrograph(scanno: int, valid=True, load_radiation_hits=True, load_audit=False, load_psf=False,
apply_rsas=False) -> xr.Dataset:
""" Load the OSIRIS Level 1 Spectrograph data for a scan.
:param int scanno: The scan number to load data for
:param bool load_radiation_hits: If True, loads radiation hit file and invalidates radiated data.
:param bool load_audit: If True, includes audit data in returned dataset
:param bool valid: If True, nans out data points that do not meet specifications
:param bool load_psf: If true, the dataset returns the wavelengths shifted by the shift parameter returned by
:class:PointSpreadFunctionFit. Also returns the fwhm for the scan.
:param bool apply_rsas: If True or 'angular', the pointing is shifted by a precomputed angular rsas value.
If 'distance', the pointing is shifted by a precomputed tangent height offset (older method).
If False (default), the pointing is not shifted.
The returned data structure contains:
<xarray.Dataset>
Dimensions: (mjd: 24, pixel: 1353)
Coordinates:
* pixel (pixel) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ...
* mjd (mjd) float64 5.239e+04 5.239e+04 5.239e+04 ...
* stw (mjd) int32 598913554 598913602 598913650 598913697 ...
Data variables:
exposureTime (mjd) float64 0.0334 0.0336 0.0336 0.0334 0.0334 ...
temperature (mjd) float64 -12.52 -12.57 -12.46 -12.54 -12.45 ...
tempavg (mjd) float64 -12.53 -12.48 -12.51 -12.47 -12.49 ...
opticstemp (mjd) float64 22.83 22.83 22.83 22.83 22.83 22.83 ...
straptemp (mjd) float64 -18.22 -18.22 -18.21 -18.21 -18.21 ...
mode (mjd) int32 4094 4094 4094 4094 4094 4094 4094 4094 ...
** The unique id code of the OSIRIS imaging mode used to collect this data. **
scienceprog (mjd) int32 21 21 21 21 21 21 21 21 21 21 21 21 21 ...
** The unique id code of the OSIRIS science mode used to collect this data. **
roe (mjd) int32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** The configuration code of the OSIRIS OS CCD readout-electronics. The pre-launch values are
defined as follows:
0 = 32x1353. No on-chip binning.
1 = 16x1353. On-chip binning of 2 reduces 32 rows to 16.
2 = 8x1353. On-chip binning of 4 reduces 32 rows to 8.
3 = 286x1353. Full CCD including storage area. Used for engineering.
4 = 143x1353. CCD imaging area. Used for engineering.
7 = 2x1353. Obsolete and unsupported mode. **
shuttermode (mjd) int32 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ...
** Indicates the mode of the OSIRIS OS shutter
0 = Close, Open, Close
1 = Closed
2 = Open **
spm_baserow (mjd) int32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** The base row used by OSIRIS in the Science Processing Module to bin the data off-chip but before
transmission to ground **
spm_numrows (mjd) int32 32 32 32 32 32 32 32 32 32 32 32 32 32 ...
** The number of rows, spatially binned in the OSIRIS Science Processing Module. **
spm_processingMode (mjd) int32 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
** The OSIRIS Science Processing Mode. The valid values are:
0 = not binned (1353xn)
1 = spatial binning (1353x1)
2 = McDade & Stegman binning (849x1)
3 = Llewellyn & Evans binning (11x32) **
targetIndex (mjd) int32 57615 57615 57615 57615 57615 57615 ...
** The real-time ACDC target index at the start of the exposure. **
exceptions (mjd) int32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** 32 bits of flags used to determine various exceptions that have occurred in processing. The most
significant bit, bit 31, indicates severity. If bit 31 is set then the record has a serious
problem and is probably unusable. If bit 31 is clear then the record has exceptions but may be
usable depending upon context. All bit fields (except bit 31) are currently t.b.d. All
definitions require consultation with level1 and level 2 OS processing groups. **
processingflags (mjd) int32 15 15 15 15 15 15 15 15 15 15 15 15 15 ...
** 32 bits of flags used to indicate which processing steps have been applied to the data.
For example data collected with the shutter closed do not have the dark current removed while all
other data do have the dark current removed. All bit fields are currently t.b.d. **
numcolumns (mjd) int32 1353 1353 1353 1353 1353 1353 1353 1353 ...
numrows (mjd) int32 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
data (mjd, pixel) float64 9.979e+10 1.286e+11 1.202e+11 ...
error (mjd, pixel) float64 3.002e+10 5.304e+10 3.555e+10 ...
flags (mjd, pixel) int32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Values of 0 mean the data is good. **
wavelength (pixel) float64 275.0 275.4 275.7 276.1 276.5 276.9 ...
** Wavelengths that match up to spectrograph pixels **
fwhm (pixel) float64 0.7763 0.7799 0.7834 0.787 0.7905 ...
** Full Width Half Max of the spectrograph point spread function **
If the audit data is loaded, the returned data structure also contains:
Coordinates:
* store (store) int64 0 1 2 3
Data variables:
audit (mjd, store) float64 5.369e+04 5.509e+04 0.0 0.0 ...
** Level 1, Audit, store=1
1.Less than 53675.0 (pre 2005-11-01)
Original software decoding.
2. 53675.0 (2005-11-01)
Fixed mjdfield to reflect the changes at Level 0.
3. 53745.0 (2006-01-10)
Applied Vega Absolute Calibration to spectrograph data. **
If the point spread function is calculated, the returned data structure also contains:
Data variables:
wavelength_error float64 0.002607
** Error in pixel / wavelength registration **
fwhm_error float64 0.005531
** Error in Full Width Half Max of the spectrograph point spread function **
For example:
>>> spectrograph = open_level1_spectrograph(scanno=6432012)
Or to get the psf adjusted wavelength and fwhm:
>>> spectrograph = open_level1_spectrograph(scanno=6432012, load_psf=True)
Or to get the raw spectrograph data:
>>> spectrograph = open_level1_spectrograph(scanno=6432012, valid=False)
"""
orbit = scanno // 1000
spectrograph = _open_level1_spectrograph_hdf4(orbit, load_audit).copy(deep=True)
spectrograph['scanno'] = scanno
# If the requested scan is 0, the previous orbit's data may need to be loaded if
# the start times between the attitude data and the spectrograph data differ.
# After merging the data from the previous orbit, trim it back to just what is needed for this orbit.
if scanno % 1000 == 0:
att = open_level1_attitude(orbit)
if spectrograph.mjd[0] - att.StartMJD.min() > MAX_ATTITUDE_DATA_TIME_DELTA:
try:
spectrograph = xr.concat(
[_open_level1_spectrograph_hdf4(orbit - 1, load_audit), spectrograph], dim='mjd')
spectrograph = spectrograph.sel(mjd=slice(att.StartMJD.min(), None))
except LookupError as e:
print('''Data for the previous orbit was expected but not found.
The attitude data begins before the spectrograph data.''')
print(e)
# Incorporate the radiation hit data
if load_radiation_hits:
radiation_hits = _load_radiation_hit_data(orbit)
radiation_hits = radiation_hits[radiation_hits.scannumber == scanno]
for _, stw, pixel, magnitude in radiation_hits.itertuples(index=False, name=None):
if magnitude < MIN_RADIATION_HIT:
spectrograph.flags[spectrograph.stw == stw, pixel] |= 0x4
else:
spectrograph.flags[spectrograph.stw == stw, pixel] |= 0x10 | 0x80
else:
warnings.warn('Radiation hit data not loaded. Radiation hits, if any, have not been marked invalid')
# Add the wavelength data and fwhm data to the pixels
spectrograph['wavelength'], spectrograph['fwhm'] = _nominal_spectrograph_wavelengths_and_fwhm()
# Replace invalid data with NaN's:
if valid:
spectrograph['data'] = spectrograph.data.where(spectrograph.flags % 2 == 0) \
.where(spectrograph.flags < 0x80) \
.where(spectrograph.exceptions == 0)
# Check for psf output
if load_psf:
from osirisl1services.psf_fit import PointSpreadFunctionFit
psf_fit = PointSpreadFunctionFit(scanno).fit()
spectrograph['fwhm'] = xr.where(
spectrograph.pixel < BLUE_RED_FWHM_PIXEL,
psf_fit.fwhm,
spectrograph['fwhm'])
spectrograph['fwhm_error'] = psf_fit.fwhm_error
spectrograph['wavelength'] = xr.where(
spectrograph.pixel < BLUE_RED_FWHM_PIXEL, # Purposefully set to the middle of the order sorter.
spectrograph['wavelength'] + psf_fit.shift,
spectrograph['wavelength'])
spectrograph['wavelength_error'] = psf_fit.shift_error
if apply_rsas == 'distance':
try:
rsas_file = Path(os.environ['ODINFLIGHTDIR'], 'rsas_meters.csv')
rsas_data = np.loadtxt(rsas_file, delimiter=',')
except OSError:
raise OSError('RSAS file ' + rsas_file + ' could not be found')
mid_day = np.array([np.nanmean(m) for m in zip(rsas_data[:, 0], rsas_data[:, 1])])
rsas_shift_m = rsas_data[:, 2]
scan_mjd_mean = float(spectrograph.mjd.mean())
# Check to make sure that there is an mjd in the rsas file close to the scan mjd
closest = mid_day[np.argmin(np.abs(mid_day - scan_mjd_mean))]
if np.abs(closest - scan_mjd_mean) > 1:
# If more than a day from the nearest rsas value raise an error
raise IndexError('Closest RSAS value in RSAS file is ' + str(np.abs(closest - scan_mjd_mean))
+ ' days away')
spectrograph['rsas_shift'] = np.interp(scan_mjd_mean, mid_day, rsas_shift_m)
spectrograph['rsas_shift'].attrs = {
'description': 'The difference in tangent altitude between the RSAS shifted and original line-of-sights',
'units': 'm'}
elif apply_rsas or apply_rsas == 'angular':
rsas = _load_angular_rasa_data()
try:
spectrograph['rsas_shift'] = rsas.sel(
mjd=spectrograph.mjd.median(), method='nearest', tolerance=1.5).item()
except KeyError:
logging.exception('RSAS, time given was not within the rsas file range.')
spectrograph['rsas_shift'].attrs = {
'description': 'The difference in angle between the RSAS shifted and original line-of-sights',
'units': 'radians'}
# Reduce the amount of data returned to the given scan
return _limit_to_scan(spectrograph, scanno)
[docs]@lru_cache(maxsize=LRU_CACHE_SIZE)
def open_level1_ir(orbit: int, channel: int, valid=True, load_audit=False) -> xr.Dataset:
""" Load the OSIRIS Level 1 IR data for an orbit and channel. Optinally load the audit.
:param int channel: IR channel number to load data from
:param int scanno: The scan number to limit the data to, default is None
:param bool valid: If True, returns the data masked to give only data points without exceptions or flags set
:param bool load_audit: If True, includes audit data in returned dataset
The returned data structure contains:
<xarray.Dataset>
Dimensions: (mjd: 43, pixel: 128)
Coordinates:
* mjd (mjd) float64 5.239e+04 5.239e+04 5.239e+04 5.239e+04 ...
* pixel (pixel) int32 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...
Data variables:
stw (mjd) int32 598913516 598913548 598913579 598913610 ...
exposureTime (mjd) float64 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ...
temperature (mjd) float64 -13.09 -12.97 -13.08 -13.21 -13.09 -13.24 ...
tempavg (mjd) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
mode (mjd) int32 1294 1294 1294 1294 1294 1294 1294 1294 ...
** The unique id code of the OSIRIS imaging mode used to collect this data. **
scienceprog (mjd) int32 4092 4092 4092 4092 4092 4092 4092 4092 ...
** The unique id code of the OSIRIS science mode used to collect this data. **
shutter (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Indicates the status of the OSIRIS IR shutter
0. Open
1. Closed **
lamp1 (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Indicates the status of the first OSIRIS IR lamp
0. Off
1. On **
lamp2 (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Indicates the status of the second OSIRIS IR lamp
0. Off
1. On **
targetIndex (mjd) int32 57615 57615 57615 57615 57615 57615 57615 ...
** The real-time ACDC target index at the start of the exposure. **
exceptions (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** 32 bits of flags used to determine various exceptions that have occurred in processing. The most
significant bit, bit 31, indicates severity. If bit 31 is set then the record has a serious
problem and is probably unusable. If bit 31 is clear then the record has exceptions but may be
usable depending upon context. All bit fields (except bit 31) are currently t.b.d. All
definitions require consultation with level1 and level 2 OS processing groups. **
processingflags (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** 32 bits of flags used to indicate which processing steps have been applied to the data.
For example data collected with the shutter closed do not have the dark current removed while all
other data do have the dark current removed. All bit fields are currently t.b.d. **
data (mjd, pixel) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
error (mjd, pixel) float64 2.3e+10 2.3e+10 2.3e+10 2.3e+10 ...
flags (mjd, pixel) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Values of 0 mean the data is good. **
If the audit data is loaded, the returned data structure also contains:
Coordinates:
* store (store) int64 0 1 2 3
Data variables:
audit (mjd, store) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
** Level 1, Audit, store=1
1.Less than 53675.0 (pre 2005-11-01)
Original software decoding.
2. 53675.0 (2005-11-01)
Fixed mjdfield to reflect the changes at Level 0.
3. 53745.0 (2006-01-10)
Applied Vega Absolute Calibration to spectrograph data. **
For example:
>>> ir = open_level1_ir(orbit=6432, channel=1)
Or if you only want the raw data:
>>> ir_raw = open_level1_ir(orbit=6432, channel=1, valid=False)
"""
assert 1 <= channel <= 3
def open_hdf4(o: int) -> xr.Dataset:
with HDF4VS(Orbit(o).l1_ir_filename) as vs:
l1 = vs[f'IR{channel} Level 1']
data = [vs['R8:R4:00128']['Array'][i[0]] for i in l1['data']]
assert_equal([vs['R8:R4:00128']['crosscheck'][i[0]] for i in l1['data']], [d[3] for d in l1['data']])
error = [vs['R8:R4:00128']['Array'][i[0]] for i in l1['error']]
assert_equal([vs['R8:R4:00128']['crosscheck'][i[0]] for i in l1['error']], [d[3] for d in l1['error']])
flags = [vs['UI1:UI1:00128']['Array'][i[0]] for i in l1['flags']]
assert_equal([vs['UI1:UI1:00128']['crosscheck'][i[0]] for i in l1['flags']], [d[3] for d in l1['flags']])
ds = xr.Dataset({'stw': (['mjd'], l1['stw']),
'exposureTime': (['mjd'], l1['exposureTime']),
'temperature': (['mjd'], l1['temperature']),
'tempavg': (['mjd'], l1['tempavg']),
# 'detectorId': (['mjd'], l1['detectorId']),
'mode': (['mjd'], l1['mode']),
'scienceprog': (['mjd'], l1['scienceprog']),
'shutter': (['mjd'], l1['shutter']),
'lamp1': (['mjd'], l1['lamp1']),
'lamp2': (['mjd'], l1['lamp2']),
'targetIndex': (['mjd'], l1['targetIndex']),
'exceptions': (['mjd'], l1['exceptions']),
'processingflags': (['mjd'], l1['processingflags']),
'data': (['mjd', 'pixel'], data),
'error': (['mjd', 'pixel'], error),
'flags': (['mjd', 'pixel'], flags)},
coords={'mjd': l1['mjd'],
'pixel': range(128)})
if load_audit:
audit = xr.Dataset({'audit': (['mjd', 'store'], l1['audit'])},
coords={'mjd': l1['mjd'],
'store': [0, 1, 2, 3]})
ds = ds.merge(audit)
return ds
ir = open_hdf4(orbit)
ir['orbit'] = orbit
ir['channel'] = channel
# If the requested scan is 0, the previous orbit's data may need to be loaded if
# the start times between the attitude data and the spectrograph data differ.
# After merging the data from the previous orbit, trim it back to just what is needed for this orbit.
att = open_level1_attitude(orbit)
if ir.mjd[0] - att.StartMJD.min() > MAX_ATTITUDE_DATA_TIME_DELTA:
try:
ir = ir.merge(open_hdf4(orbit - 1))
ir = ir.sel(mjd=slice(att.StartMJD.min(), None))
except LookupError as e:
print('''Data for the previous orbit was expected but not found.
The attitude data begins before the spectrograph data.''')
print(e)
# Mask ir to return valid data points if valid is set to true
if valid:
ir['data'] = ir.data.where(ir.exceptions == 0).where(ir.flags == 0)
return ir
@lru_cache(maxsize=LRU_CACHE_SIZE)
def open_level1b_ir(orbit: Orbit, channel: int, valid=True, load_audit=False) -> xr.Dataset:
""" Load the OSIRIS Level 1b IR data for an orbit and channel. Optinally load the audit.
The returned data structure contains:
<xarray.Dataset>
Dimensions: (mjd: 30, pixel: 128)
Coordinates:
* mjd (mjd) float64 5.239e+04 5.239e+04 5.239e+04 5.239e+04 ...
* pixel (pixel) int32 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...
Data variables:
stw (mjd) int32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
exposureTime (mjd) float64 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ...
temperature (mjd) float64 -13.09 -12.97 -13.08 -13.21 -13.09 -13.24 ...
tempavg (mjd) float64 -9.999e+03 -9.999e+03 -9.999e+03 ...
detectorId (mjd) int32 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
mode (mjd) int32 1294 1294 1294 1294 1294 1294 1294 1294 ...
** The unique id code of the OSIRIS imaging mode used to collect this data. **
scienceprog (mjd) int32 4092 4092 4092 4092 4092 4092 4092 4092 ...
** The unique id code of the OSIRIS science mode used to collect this data. **
shutter (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Indicates the status of the OSIRIS IR shutter
0. Open
1. Closed **
lamp1 (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Indicates the status of the first OSIRIS IR lamp
0. Off
1. On **
lamp2 (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** Indicates the status of the second OSIRIS IR lamp
0. Off
1. On **
targetIndex (mjd) int32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** The real-time ACDC target index at the start of the exposure. **
exceptions (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** 32 bits of flags used to determine various exceptions that have occurred in processing. The most
significant bit, bit 31, indicates severity. If bit 31 is set then the record has a serious
problem and is probably unusable. If bit 31 is clear then the record has exceptions but may be
usable depending upon context. All bit fields (except bit 31) are currently t.b.d. All
definitions require consultation with level1 and level 2 OS processing groups. **
processingflags (mjd) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
** 32 bits of flags used to indicate which processing steps have been applied to the data.
For example data collected with the shutter closed do not have the dark current removed while all
other data do have the dark current removed. All bit fields are currently t.b.d. **
data (mjd, pixel) float64 -9.999e+03 -9.999e+03 -9.999e+03 ...
error (mjd, pixel) float64 -9.999e+03 -9.999e+03 -9.999e+03 ...
flags (mjd, pixel) int32 128 128 128 128 128 128 128 128 128 ...
** Values of 0 mean the data is good. **
If the audit data is loaded, the returned data structure also contains:
Coordinates:
* store (store) int64 0 1 2 3
Data variables:
audit (mjd, store) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
** Level 1, Audit, store=1
1.Less than 53675.0 (pre 2005-11-01)
Original software decoding.
2. 53675.0 (2005-11-01)
Fixed mjdfield to reflect the changes at Level 0.
3. 53745.0 (2006-01-10)
Applied Vega Absolute Calibration to spectrograph data. **
For example:
>>> irb = open_level1b_ir(orbit=6432, channel=1)
Or, if you only want the raw data:
>>> irb_raw = open_level1b_ir(orbit=6432, channel=1, valid=False)
"""
assert channel == 1
def open_hdf4(o: int) -> xr.Dataset:
with HDF4VS(Orbit(o).l1b_ir_filename) as vs:
l1 = vs[f'IR{channel} Level 1']
data = [vs['R8:R4:00128']['Array'][i[0]] for i in l1['data']]
assert_equal([vs['R8:R4:00128']['crosscheck'][i[0]] for i in l1['data']], [d[3] for d in l1['data']])
error = [vs['R8:R4:00128']['Array'][i[0]] for i in l1['error']]
assert_equal([vs['R8:R4:00128']['crosscheck'][i[0]] for i in l1['error']], [d[3] for d in l1['error']])
flags = [vs['UI1:UI1:00128']['Array'][i[0]] for i in l1['flags']]
assert_equal([vs['UI1:UI1:00128']['crosscheck'][i[0]] for i in l1['flags']], [d[3] for d in l1['flags']])
ds = xr.Dataset({'stw': (['mjd'], l1['stw']),
'exposureTime': (['mjd'], l1['exposureTime']),
'temperature': (['mjd'], l1['temperature']),
'tempavg': (['mjd'], l1['tempavg']),
'detectorId': (['mjd'], l1['detectorId']),
'mode': (['mjd'], l1['mode']),
'scienceprog': (['mjd'], l1['scienceprog']),
'shutter': (['mjd'], l1['shutter']),
'lamp1': (['mjd'], l1['lamp1']),
'lamp2': (['mjd'], l1['lamp2']),
'targetIndex': (['mjd'], l1['targetIndex']),
'exceptions': (['mjd'], l1['exceptions']),
'processingflags': (['mjd'], l1['processingflags']),
'data': (['mjd', 'pixel'], data),
'error': (['mjd', 'pixel'], error),
'flags': (['mjd', 'pixel'], flags)},
coords={'mjd': l1['mjd'],
'pixel': range(128)})
if load_audit:
audit = xr.Dataset({'audit': (['mjd', 'store'], l1['audit'])},
coords={'mjd': l1['mjd'],
'store': [0, 1, 2, 3]})
ds = ds.merge(audit)
return ds
ir = open_hdf4(orbit)
ir['orbit'] = orbit
ir['channel'] = channel
# If the whole orbit is requested or the requested scan is 0, the previous orbit's data may need to be loaded if
# the start times between the attitude data and the spectrograph data differ.
# After merging the data from the previous orbit, trim it back to just what is needed for this orbit.
att = open_level1_attitude(orbit)
if ir.mjd[0] - att.StartMJD.min() > MAX_ATTITUDE_DATA_TIME_DELTA:
try:
ir = ir.merge(open_hdf4(orbit - 1))
ir = ir.sel(mjd=slice(att.StartMJD.min(), None))
except LookupError as e:
print('''Data for the previous orbit was expected but not found.
The attitude data begins before the spectrograph data.''')
print(e)
# Return only data without exceptions or flags by applying a mask
if valid:
ir['data'] = ir.data.where(ir.exceptions == 0).where(ir.flags == 0)
return ir
@lru_cache(maxsize=LRU_CACHE_SIZE)
def open_level1_ecmwf_orbit(orbit: int) -> xr.Dataset:
""" Load the OSIRIS Level 1 ECMWF data for an orbit.
:param int orbit: Orbit number for the data
The returned data structure contains:
<xarray.Dataset>
Dimensions: (Altitude: 201, ScanNumber: 66)
Coordinates:
* ScanNumber (ScanNumber) int64 6432000 6432001 6432002 ... 6432064 6432065
* Altitude (Altitude) float64 0.0 500.0 1e+03 ... 9.9e+04 9.95e+04 1e+05
Data variables:
density (ScanNumber, Altitude) float64 2.469e+19 ... 1.021e+13
** Air Density **
temperature (ScanNumber, Altitude) float64 261.1 260.3 261.4 261.7 ...
** Air Temperature **
For example:
>>> osiris_ecmwf = open_level1_ecmwf(6432)
"""
with HDF4VS(Orbit(orbit).l1_ecmwf_filename) as vs:
ecmwf = vs['Osiris ECMWF']
assert (all(alt == 0.0 for alt in ecmwf['ecmwf_minalt']))
assert (all(alt == 100000.0 for alt in ecmwf['ecmwf_maxalt']))
# Entries of [0, 0, 0, 0] in the data pointers represent no data
# Disregard those scans from the asserts and drop correspondng scans from the returned dataset
altitude = [vs['R8:R8:00221']['Array'][i[0]] for i in ecmwf['altitudes']]
assert_equal([vs['R8:R8:00221']['crosscheck'][i[0]] for i in ecmwf['altitudes'] if i != [0, 0, 0, 0]],
[d[3] for d in ecmwf['altitudes'] if d != [0, 0, 0, 0]])
# Ensure all the altitudes range between 0 km and 110 km at 500m intervals
for column in altitude:
assert_array_equal(column, np.arange(0.0, 110500.0, 500.0))
density = [vs['R8:R8:00221']['Array'][i[0]] for i in ecmwf['Density']]
assert_equal([vs['R8:R8:00221']['crosscheck'][i[0]] for i in ecmwf['Density'] if i != [0, 0, 0, 0]],
[d[3] for d in ecmwf['Density'] if d != [0, 0, 0, 0]])
temperature = [vs['R8:R8:00221']['Array'][i[0]] for i in ecmwf['Temperature']]
assert_equal([vs['R8:R8:00221']['crosscheck'][i[0]] for i in ecmwf['Temperature'] if i != [0, 0, 0, 0]],
[d[3] for d in ecmwf['Temperature'] if d != [0, 0, 0, 0]])
# Check that all of the no data [0, 0, 0, 0] entries were in the same positions
for (a, d, t) in zip(ecmwf['altitudes'], ecmwf['Density'], ecmwf['Temperature']):
if a == [0, 0, 0, 0]:
assert (d == [0, 0, 0, 0]) and (t == [0, 0, 0, 0])
ds = xr.Dataset({'density': (['ScanNumber', 'Altitude'], density),
'temperature': (['ScanNumber', 'Altitude'], temperature)},
coords={'ScanNumber': ecmwf['ScanNumber'],
'Altitude': altitude[0]})
# Drop scans with no data and return only up to 100km
ds = ds.sel(ScanNumber=[i != [0, 0, 0, 0] for i in ecmwf['altitudes']], Altitude=slice(0.0, 100000.0))
return ds
[docs]def open_level1_ecmwf(scanno: int) -> xr.Dataset:
""" Load the OSIRIS Level 1 ECMWF data for a scan.
:param int scanno: Scan number for the data
The returned data structure contains:
<xarray.Dataset>
Dimensions: (Altitude: 201)
Coordinates:
* ScanNumber int32 6432012
* Altitude (Altitude) float64 0.0 500.0 1e+03 1.5e+03 2e+03 2.5e+03 ...
Data variables:
density (ScanNumber, Altitude) float64 2.839e+19 2.652e+19 .....
** Air Density **
temperature (ScanNumber, Altitude) float64 261.1 260.3 261.4 261.7 ...
** Air Temperature **
For example:
>>> osiris_ecmwf = open_level1_ecmwf(6432012)
"""
orbit = scanno // 1000
ds = open_level1_ecmwf_orbit(orbit)
return ds.sel(ScanNumber=scanno)