From c44da0bd20c9083f4361a315f5f7e0b6c7a1945d Mon Sep 17 00:00:00 2001 From: Jesse Mapel Date: Thu, 25 Jul 2019 10:40:52 -0700 Subject: [PATCH] updates data_isis to return positions in the correct frame (#214) * Added helper func to convert ISIS table to rotations * Added ISIS frame chain * Fixed rotation direction in data_isis * First pass at data_isis sun position * added sensor position to data_isis * More updates * Updated data_isis test after removing label * Updated data_isis to new frame chain * updated to get tests passing * Cleaned up how ISIS tables are read * Fixed missing return * Cleaned up fixtures in test_transformation * Added km to m conversion --- ale/base/data_isis.py | 346 +++++++++++---------------- ale/encoders.py | 2 + ale/rotation.py | 16 +- ale/transformation.py | 79 ++++++ tests/pytests/test_data_isis.py | 123 +++++++++- tests/pytests/test_rotation.py | 9 + tests/pytests/test_transformation.py | 56 ++++- 7 files changed, 415 insertions(+), 216 deletions(-) diff --git a/ale/base/data_isis.py b/ale/base/data_isis.py index e4808c5..3fb6b2f 100644 --- a/ale/base/data_isis.py +++ b/ale/base/data_isis.py @@ -4,15 +4,16 @@ import struct from glob import glob import numpy as np +from numpy.polynomial.polynomial import polyval, polyder from dateutil import parser import pvl import spiceypy as spice +from ale.rotation import ConstantRotation, TimeDependentRotation +from ale.transformation import FrameChain from ale import config -from scipy.interpolate import interp1d -from scipy.spatial.transform import Slerp -from scipy.spatial.transform import Rotation +from scipy.interpolate import interp1d, BPoly from ale.base.label_isis import IsisLabel @@ -36,187 +37,145 @@ def read_table_data(table_label, cube): cubehandle.seek(table_label['StartByte'] - 1) return cubehandle.read(table_label['Bytes']) -def field_size(field_label): +def parse_table(table_label, data): """ - Helper function to determine the size of a binary - table field - - Parameters - ---------- - field_label : PVLModule - The field label - - Returns - ------- - int : - The size of the one entry in bytes - """ - data_sizes = { - 'Integer' : 4, - 'Double' : 8, - 'Real' : 4, - 'Text' : 1 - } - return data_sizes[field_label['Type']] * field_label['Size'] - -def field_format(field_label): - """ - Helper function to get the format string for a - single entry in a table field - - Parameters - ---------- - field_label : PVLModule - The field label - - Returns - ------- - str : - The format string for the entry binary - """ - data_formats = { - 'Integer' : 'i', - 'Double' : 'd', - 'Real' : 'f' - } - return data_formats[field_label['Type']] * field_label['Size'] - -def parse_field(field_label, data, encoding='latin_1'): - """ - Parses a binary table field entry and converts it into - an in memory data type - - Parameters - ---------- - field_label : PVLModule - The field label - - data : bytes - The binary data for the field entry - - Returns - ------- - Union[int, float, str, list] : - The table field entry converted to a native python - type - """ - if field_label['Type'] == 'Text': - field_data = data[:field_label['Size']].decode(encoding=encoding) - else: - data_format = field_format(field_label) - field_data = struct.unpack_from(data_format, data) - if len(field_data) == 1: - field_data = field_data[0] - return field_data - -def parse_table_data(table_label, data): - """ - Parses an ISIS table into a dict where the keys are the - field names and the values are lists of entries. + Parse an ISIS table into a dictionary. Parameters ---------- table_label : PVLModule - The table label - + The ISIS table label data : bytes - The binary data for the entire table + The binary component of the ISIS table Returns ------- dict : - The table as a dict + The table as a dictionary with the keywords from the label and the + binary data """ + data_sizes = {'Integer' : 4, + 'Double' : 8, + 'Real' : 4, + 'Text' : 1} + data_formats = {'Integer' : 'i', + 'Double' : 'd', + 'Real' : 'f'} + + # Parse the binary data fields = table_label.getlist('Field') results = {field['Name']:[] for field in fields} offset = 0 for record in range(table_label['Records']): for field in fields: - field_data = parse_field(field, data[offset:]) + if field['Type'] == 'Text': + field_data = data[offset:offset+field['Size']].decode(encoding='latin_1') + else: + data_format = data_formats[field['Type']] * field['Size'] + field_data = struct.unpack_from(data_format, data[offset:]) + if len(field_data) == 1: + field_data = field_data[0] + results[field['Name']].append(field_data) - offset += field_size(field) + offset += data_sizes[field['Type']] * field['Size'] + + # Parse the keywords from the label + results.update({key : value for key, value in table_label.items() if not isinstance(value, pvl._collections.PVLGroup)}) + return results -def parse_rotation_table(label, field_data): +def rotate_positions(table, rotation): """ - Parses ISIS rotation table data. + Rotate the positions in a position Table. + + If the table stores positions as a function, then it will re compute them + based on the original size of the table. Parameters ---------- - label : PVLModule - The table label - - field_data : dict - The table data as a dict with field names - as keys and lists of entries as values + table : dict + The position table as a dictionary + rotation : TimeDependentRotation + The rotation to rotate the positions by Returns ------- - dict : - The rotation data - """ - results = {} - if all (key in field_data for key in ('J2000Q0','J2000Q1','J2000Q2','J2000Q3')): - results['Rotations'] = [ [q0, q1, q2, q3] for q0, q1, q2, q3 in zip(field_data['J2000Q0'],field_data['J2000Q1'],field_data['J2000Q2'],field_data['J2000Q3']) ] - if all (key in field_data for key in ('AV1','AV2','AV3')): - results['AngularVelocities'] = np.array( [ [av1, av2, av3] for av1, av2, av3 in zip(field_data['AV1'],field_data['AV2'],field_data['AV3']) ] ) - if 'ET' in field_data: - results['Times'] = np.array(field_data['ET']) - if all (key in field_data for key in ('J2000Ang1','J2000Ang2','J2000Ang3')): - results['EulerCoefficients'] = np.array([field_data['J2000Ang1'],field_data['J2000Ang2'],field_data['J2000Ang3']]) - results['BaseTime'] = field_data['J2000Ang1'][-1] - results['TimeScale'] = field_data['J2000Ang2'][-1] - - if 'TimeDependentFrames' in label: - results['TimeDependentFrames'] = np.array(label['TimeDependentFrames']) - if 'CkTableOriginalSize' in label: - results['CkTableOriginalSize'] = label['CkTableOriginalSize'] - if all (key in label for key in ('ConstantRotation','ConstantFrames')): - const_rotation_mat = np.array(label['ConstantRotation']) - results['ConstantRotation'] = np.reshape(const_rotation_mat, (3, 3)) - results['ConstantFrames'] = np.array(label['ConstantFrames']) - if all (key in label for key in ('PoleRa','PoleDec','PrimeMeridian')): - results['BodyRotationCoefficients'] = np.array( [label['PoleRa'],label['PoleDec'],label['PrimeMeridian']] ) - if all (key in label for key in ('PoleRaNutPrec','PoleDecNutPrec','PmNutPrec','SysNutPrec0','SysNutPrec1')): - results['SatelliteNutationPrecessionCoefficients'] = np.array( [label['PoleRaNutPrec'],label['PoleDecNutPrec'],label['PmNutPrec']] ) - results['PlanetNutationPrecessionAngleCoefficients'] = np.array( [label['SysNutPrec0'],label['SysNutPrec1']] ) + : 2darray + Array of rotated positions + : array + Array of times for the positions - return results + """ + # Case 1, the table has positions at discrete times + if 'J2000X' in table: + ephemeris_times = table['ET'] + positions = 1000 * np.array([table['J2000X'], + table['J2000Y'], + table['J2000Z']]).T + rotated_pos = rotation.apply_at(positions, ephemeris_times) + # Case 2, the table has coefficients of polynomials for the position + elif 'J2000SVX' in table: + ephemeris_times = np.linspace(table['SpkTableStartTime'], + table['SpkTableEndTime'], + table['SpkTableOriginalSize']) + base_time = table['J2000SVX'][-1] + time_scale = table['J2000SVY'][-1] + scaled_times = (ephemeris_times - base_time) / time_scale + coeffs = np.array([table['J2000SVX'][:-1], + table['J2000SVY'][:-1], + table['J2000SVZ'][:-1]]).T + positions = 1000 * polyval(scaled_times, coeffs).T + rotated_pos = rotation.apply_at(positions, ephemeris_times) + else: + raise ValueError('No positions are available in the input table.') + return rotated_pos, ephemeris_times -def parse_position_table(label, field_data): +def rotate_velocities(table, rotation): """ - Parses ISIS position table data. + Rotate the velocities in a position Table. + + If the table stores positions as a function, then it will re compute them + based on the original size of the table. Parameters ---------- - label : PVLModule - The table label - - field_data : dict - The table data as a dict with field names as keys - and lists of entries as values + table : dict + The position table as a dict from parse_position_table + rotation : TimeDependentRotation + The rotation to rotate the velocities by Returns ------- - dict : - The position data - """ - results = {} - if all (key in field_data for key in ('J2000X','J2000Y','J2000Z')): - results['Positions'] = np.array( [ [x, y, z] for x, y, z in zip(field_data['J2000X'],field_data['J2000Y'],field_data['J2000Z']) ] ) - if 'ET' in field_data: - results['Times'] = np.array(field_data['ET']) - if all (key in field_data for key in ('J2000XV','J2000YV','J2000ZV')): - results['Velocities'] = np.array( [ [x, y, z] for x, y, z in zip(field_data['J2000XV'],field_data['J2000YV'],field_data['J2000ZV']) ] ) - if all (key in field_data for key in ('J2000SVX','J2000SVY','J2000SVZ')): - results['PositionCoefficients'] = np.array( [field_data['J2000SVX'][:-1],field_data['J2000SVY'][:-1],field_data['J2000SVZ'][:-1]] ) - results['BaseTime'] = field_data['J2000SVX'][-1] - results['TimeScale'] = field_data['J2000SVY'][-1] - - if 'SpkTableOriginalSize' in label: - results['SpkTableOriginalSize'] = label['SpkTableOriginalSize'] - return results + : 2darray + Array of rotated velocities. Returns None if no velocities are in the table. + """ + # Case 1, the table has velocities at discrete times + if 'J2000XV' in table: + ephemeris_times = table['ET'] + velocity = 1000 * np.array([table['J2000XV'], + table['J2000YV'], + table['J2000ZV']]).T + rotated_vel = rotation.apply_at(velocity, ephemeris_times) + # Case 2, the table has coefficients of polynomials for the position + elif 'J2000SVX' in table: + ephemeris_times = np.linspace(table['SpkTableStartTime'], + table['SpkTableEndTime'], + table['SpkTableOriginalSize']) + base_time = table['J2000SVX'][-1] + time_scale = table['J2000SVY'][-1] + scaled_times = (ephemeris_times - base_time) / time_scale + coeffs = np.array([table['J2000SVX'][:-1], + table['J2000SVY'][:-1], + table['J2000SVZ'][:-1]]) + scaled_vel = 1000 * polyval(scaled_times, polyder(coeffs,axis=1).T).T + # We took a derivative in scaled time, so we have to multiply by our + # scale in order to get the derivative in real time + velocity = scaled_vel / time_scale + rotated_vel = rotation.apply_at(velocity, ephemeris_times) + else: + rotated_vel = None + return rotated_vel class IsisSpice(): """Mixin class for reading from an ISIS cube that has been spiceinit'd @@ -259,24 +218,6 @@ class IsisSpice(): """ - @property - def label(self): - """ - Loads a PVL from from the _file attribute and - parses the binary table data. - - Returns - ------- - PVLModule : - Dict-like object with PVL keys - """ - if not hasattr(self, "_label"): - try: - self._label = pvl.load(self.file) - except: - raise ValueError("{} is not a valid label".format(self.file)) - return self._label - @property def inst_pointing_table(self): """ @@ -292,8 +233,7 @@ class IsisSpice(): for table in self.label.getlist('Table'): if table['Name'] == 'InstrumentPointing': binary_data = read_table_data(table, self._file) - field_data = parse_table_data(table, binary_data) - self._inst_pointing_table = parse_rotation_table(table, field_data) + self._inst_pointing_table = parse_table(table, binary_data) return self._inst_pointing_table @property @@ -311,8 +251,7 @@ class IsisSpice(): for table in self.label.getlist('Table'): if table['Name'] == 'BodyRotation': binary_data = read_table_data(table, self._file) - field_data = parse_table_data(table, binary_data) - self._body_orientation_table = parse_rotation_table(table, field_data) + self._body_orientation_table = parse_table(table, binary_data) return self._body_orientation_table @property @@ -330,8 +269,7 @@ class IsisSpice(): for table in self.label.getlist('Table'): if table['Name'] == 'InstrumentPosition': binary_data = read_table_data(table, self._file) - field_data = parse_table_data(table, binary_data) - self._inst_position_table = parse_position_table(table, field_data) + self._inst_position_table = parse_table(table, binary_data) return self._inst_position_table @property @@ -349,8 +287,7 @@ class IsisSpice(): for table in self.label.getlist('Table'): if table['Name'] == 'SunPosition': binary_data = read_table_data(table, self._file) - field_data = parse_table_data(table, binary_data) - self._sun_position_table = parse_position_table(table, field_data) + self._sun_position_table = parse_table(table, binary_data) return self._sun_position_table def __enter__(self): @@ -571,6 +508,25 @@ class IsisSpice(): if re.match(regex, key[0]): return self.isis_naif_keywords[key[0]] + @property + def frame_chain(self): + """ + Return the root node of the rotation frame tree/chain. + + The root node is the J2000 reference frame. The other nodes in the + tree can be accessed via the methods in the FrameNode class. + + Returns + ------- + FrameNode + The root node of the frame tree. This will always be the J2000 reference frame. + """ + if not hasattr(self, '_frame_chain'): + self._frame_chain = FrameChain.from_isis_tables( + inst_pointing = self.inst_pointing_table, + body_orientation = self.body_orientation_table) + return self._frame_chain + @property def sun_position(self): @@ -587,9 +543,10 @@ class IsisSpice(): of the target body in the J2000 reference frame as a tuple of numpy arrays. """ - return (self.sun_position_table.get('Positions', 'None'), - self.sun_position_table.get('Velocities', 'None'), - self.sun_position_table.get('Times', 'None')) + j2000_to_target = self.frame_chain.compute_rotation(1, self.target_frame_id) + positions, times = rotate_positions(self.sun_position_table, j2000_to_target) + velocities = rotate_velocities(self.sun_position_table, j2000_to_target) + return positions, velocities, times @property @@ -607,30 +564,10 @@ class IsisSpice(): : (positions, velocities, times) a tuple containing a list of positions, a list of velocities, and a list of times """ - inst_positions_times = np.linspace(self.inst_position_table["Times"][0], - self.inst_position_table["Times"][-1], - self.number_of_ephemerides) - - # interpolate out positions - if len(self.inst_position_table["Times"]) < 2: - time_0 = self.inst_position_table["Times"][0] - position_0 = self.inst_position_table["Positions"][0] - velocity_0 = self.inst_position_table["Velocities"][0] - coefs = np.asarray([position_0 - time_0*velocity_0, - velocity_0]) - positions = np.polynomial.polynomial.polyval(inst_positions_times, coefs) - - else: - f_positions_x = interp1d(self.inst_position_table["Times"], self.inst_position_table["Positions"].T[0]) - f_positions_y = interp1d(self.inst_position_table["Times"], self.inst_position_table["Positions"].T[1]) - f_positions_z = interp1d(self.inst_position_table["Times"], self.inst_position_table["Positions"].T[2]) - - positions = np.asarray([f_positions_x(inst_positions_times), - f_positions_y(inst_positions_times), - f_positions_z(inst_positions_times)]) - - # convert positions to Body-Fixed and scale to meters - return self._body_j2k2bf_rotation.apply(positions.T)*1000 + j2000_to_target = self.frame_chain.compute_rotation(1, self.target_frame_id) + positions, times = rotate_positions(self.inst_position_table, j2000_to_target) + velocities = rotate_velocities(self.inst_position_table, j2000_to_target) + return positions, velocities, times @property @@ -663,7 +600,6 @@ class IsisSpice(): """ return self.isis_naif_keywords["INS{}_OD_K".format(self.ikid)] - @property def sensor_frame_id(self): if 'ConstantFrames' in self.inst_pointing_table: @@ -674,7 +610,7 @@ class IsisSpice(): @property def target_frame_id(self): - if 'ConstantFrames' in self.inst_pointing_table: - return self.body_rotation_table['ConstantFrames'][0] + if 'ConstantFrames' in self.body_orientation_table: + return self.body_orientation_table['ConstantFrames'][0] else: - return self.body_rotation_table['TimeDependentFrames'][0] + return self.body_orientation_table['TimeDependentFrames'][0] diff --git a/ale/encoders.py b/ale/encoders.py index a7e5c1f..d30c62e 100644 --- a/ale/encoders.py +++ b/ale/encoders.py @@ -4,6 +4,8 @@ import datetime class NumpyEncoder(json.JSONEncoder): def default(self, obj): + if isinstance(obj, set): + return list(obj) if isinstance(obj, np.integer): return int(obj) elif isinstance(obj, np.floating): diff --git a/ale/rotation.py b/ale/rotation.py index b17bd34..8b8f6b0 100644 --- a/ale/rotation.py +++ b/ale/rotation.py @@ -127,7 +127,7 @@ class TimeDependentRotation: The NAIF ID code for the destination frame """ - def from_euler(sequence, euler, times, source, dest): + def from_euler(sequence, euler, times, source, dest, degrees=False): """ Create a time dependent rotation from a set of Euler angles. @@ -145,12 +145,15 @@ class TimeDependentRotation: The NAIF ID code for the source frame dest : int The NAIF ID code for the destination frame + degrees : bool + If the angles are in degrees. If false, then degrees are + assumed to be in radians. Defaults to False. See Also -------- scipy.spatial.transform.Rotation.from_euler """ - rot = Rotation.from_euler(sequence, np.asarray(euler)) + rot = Rotation.from_euler(sequence, np.asarray(euler), degrees=degrees) return TimeDependentRotation(rot.as_quat(), times, source, dest) def __init__(self, quats, times, source, dest): @@ -245,3 +248,12 @@ class TimeDependentRotation: return TimeDependentRotation(new_quats, new_times, other.source, self.dest) else: raise TypeError("Rotations can only be composed with other rotations.") + + def apply_at(self, vec, et): + """ + Apply the rotation at a specific time + """ + if len(self.times) == 1 and self.times[0] == et: + return self._rots.apply(vec) + rot = Slerp(self.times, self._rots)(et) + return rot.apply(vec) diff --git a/ale/transformation.py b/ale/transformation.py index 1a9263e..88bb6de 100644 --- a/ale/transformation.py +++ b/ale/transformation.py @@ -1,4 +1,5 @@ import numpy as np +from numpy.polynomial.polynomial import polyval, polyder import networkx as nx from networkx.algorithms.shortest_paths.generic import shortest_path @@ -6,6 +7,69 @@ import spiceypy as spice from ale.rotation import ConstantRotation, TimeDependentRotation +def create_rotations(rotation_table): + """ + Convert an ISIS rotation table into rotation objects. + + Parameters + ---------- + rotation_table : dict + The rotation ISIS table as a dictionary + + Returns + ------- + : list + A list of time dependent or constant rotation objects from the table. This + list will always have either 1 or 2 elements. The first rotation will be + time dependent and the second rotation will be constant. The rotations will + be ordered such that the reference frame the first rotation rotates to is + the reference frame the second rotation rotates from. + """ + rotations = [] + root_frame = rotation_table['TimeDependentFrames'][-1] + last_time_dep_frame = rotation_table['TimeDependentFrames'][0] + # Case 1: It's a table of quaternions and times + if 'J2000Q0' in rotation_table: + # SPICE quaternions are (W, X, Y, Z) and ALE uses (X, Y, Z, W). + quats = np.array([rotation_table['J2000Q1'], + rotation_table['J2000Q2'], + rotation_table['J2000Q3'], + rotation_table['J2000Q0']]).T + time_dep_rot = TimeDependentRotation(quats, + rotation_table['ET'], + root_frame, + last_time_dep_frame) + rotations.append(time_dep_rot) + # Case 2: It's a table of Euler angle coefficients + elif 'J2000Ang1' in rotation_table: + ephemeris_times = np.linspace(rotation_table['CkTableStartTime'], + rotation_table['CkTableEndTime'], + rotation_table['CkTableOriginalSize']) + base_time = rotation_table['J2000Ang1'][-1] + time_scale = rotation_table['J2000Ang2'][-1] + scaled_times = (ephemeris_times - base_time) / time_scale + coeffs = np.array([rotation_table['J2000Ang1'][:-1], + rotation_table['J2000Ang2'][:-1], + rotation_table['J2000Ang3'][:-1]]).T + angles = polyval(scaled_times, coeffs).T + # ISIS is hard coded to ZXZ (313) Euler angle axis order. + time_dep_rot = TimeDependentRotation.from_euler('zxz', + angles, + ephemeris_times, + root_frame, + last_time_dep_frame, + degrees=True) + rotations.append(time_dep_rot) + + if 'ConstantRotation' in rotation_table: + last_constant_frame = rotation_table['ConstantFrames'][0] + rot_mat = np.reshape(np.array(rotation_table['ConstantRotation']), (3, 3)) + constant_rot = ConstantRotation.from_matrix(rot_mat, + last_time_dep_frame, + last_constant_frame) + rotations.append(constant_rot) + return rotations + class FrameChain(nx.DiGraph): """ This class is responsible for handling rotations between reference frames. @@ -40,6 +104,21 @@ class FrameChain(nx.DiGraph): frame_chain.add_edge(s, d, rotation=rotation) return frame_chain + @classmethod + def from_isis_tables(cls, *args, inst_pointing = {}, body_orientation={}, **kwargs): + frame_chain = cls() + + for rotation in create_rotations(inst_pointing): + frame_chain.add_edge(rotation.source, + rotation.dest, + rotation=rotation) + + for rotation in create_rotations(body_orientation): + frame_chain.add_edge(rotation.source, + rotation.dest, + rotation=rotation) + return frame_chain + def add_edge(self, s, d, rotation, **kwargs): super(FrameChain, self).add_edge(s, d, rotation=rotation, **kwargs) super(FrameChain, self).add_edge(d, s, rotation=rotation.inverse(), **kwargs) diff --git a/tests/pytests/test_data_isis.py b/tests/pytests/test_data_isis.py index 10754e7..ad66634 100644 --- a/tests/pytests/test_data_isis.py +++ b/tests/pytests/test_data_isis.py @@ -1,5 +1,6 @@ import pytest import pvl +import numpy as np from ale.base.data_isis import IsisSpice @@ -503,7 +504,7 @@ End @pytest.fixture def testdata(): isis_spice = IsisSpice() - isis_spice._label = pvl.loads(testlabel) + isis_spice.label = pvl.loads(testlabel) return isis_spice @@ -570,3 +571,123 @@ def test_isis_naif_keywords(testdata): def test_odtk(testdata): assert testdata.odtk == [-0.0048509, 2.41312e-07, -1.62369e-13] + + +def test_frame_chain(testdata): + testdata._inst_pointing_table = { + 'J2000Q0' : [1, 0.5], + 'J2000Q1' : [0, 0.5], + 'J2000Q2' : [0, 0.5], + 'J2000Q3' : [0, 0.5], + 'ET' : [0, 1], + 'TimeDependentFrames' : [-1000, -100, 1], + 'ConstantRotation' : [0, 0, 1, 1, 0 , 0, 0, 1, 0], + 'ConstantFrames' : [-1020, -1000]} + testdata._body_orientation_table = { + 'J2000Q0' : [1, 0.5], + 'J2000Q1' : [0, 0.5], + 'J2000Q2' : [0, 0.5], + 'J2000Q3' : [0, 0.5], + 'ET' : [0, 1], + 'TimeDependentFrames' : [80, 1], + 'ConstantRotation' : [0, 0, 1, 1, 0 , 0, 0, 1, 0], + 'ConstantFrames' : [81, 80]} + frame_chain = testdata.frame_chain + assert len(frame_chain.nodes) == 5 + assert frame_chain.has_node(1) + assert frame_chain.has_node(80) + assert frame_chain.has_node(81) + assert frame_chain.has_node(-1000) + assert frame_chain.has_node(-1020) + assert len(frame_chain.edges) == 8 + assert frame_chain.has_edge(1, 80) + assert frame_chain.has_edge(80, 1) + assert frame_chain.has_edge(81, 80) + assert frame_chain.has_edge(80, 81) + assert frame_chain.has_edge(1, -1000) + assert frame_chain.has_edge(-1000, 1) + assert frame_chain.has_edge(-1020, -1000) + assert frame_chain.has_edge(-1000, -1020) + +def test_sun_position_cache(testdata): + testdata._inst_pointing_table = { + 'J2000Q0' : [1, 1], + 'J2000Q1' : [0, 0], + 'J2000Q2' : [0, 0], + 'J2000Q3' : [0, 0], + 'ET' : [0, 1], + 'TimeDependentFrames' : [-1000, -100, 1]} + testdata._body_orientation_table = { + 'J2000Q0' : [1, 0.5], + 'J2000Q1' : [0, 0.5], + 'J2000Q2' : [0, 0.5], + 'J2000Q3' : [0, 0.5], + 'ET' : [0, 1], + 'TimeDependentFrames' : [80, 1]} + testdata._sun_position_table = { + 'J2000X' : [1, 0], + 'J2000Y' : [0, 1], + 'J2000Z' : [0, 0], + 'J2000XV' : [-1, 0], + 'J2000YV' : [0, -1], + 'J2000ZV' : [0, 0], + 'ET' : [0, 1]} + sun_pos, sun_vel, sun_times = testdata.sun_position + np.testing.assert_almost_equal(sun_pos, np.array([[1000, 0, 0], [0, 0, 1000]])) + np.testing.assert_almost_equal(sun_vel, np.array([[-1000, 0, 0], [0, 0, -1000]])) + np.testing.assert_equal(sun_times, np.array([0, 1])) + +def test_sun_position_polynomial(testdata): + testdata._inst_pointing_table = { + 'J2000Q0' : [1, 1], + 'J2000Q1' : [0, 0], + 'J2000Q2' : [0, 0], + 'J2000Q3' : [0, 0], + 'ET' : [2, 4], + 'TimeDependentFrames' : [-1000, -100, 1]} + testdata._body_orientation_table = { + 'J2000Q0' : [1, 0.5], + 'J2000Q1' : [0, 0.5], + 'J2000Q2' : [0, 0.5], + 'J2000Q3' : [0, 0.5], + 'ET' : [2, 4], + 'TimeDependentFrames' : [80, 1]} + testdata._sun_position_table = { + 'SpkTableOriginalSize' : 2, + 'SpkTableStartTime' : 2, + 'SpkTableEndTime' : 4, + 'J2000SVX' : [1, -1, 2], + 'J2000SVY' : [0, 1, 2], + 'J2000SVZ' : [0, -1, 1]} + sun_pos, sun_vel, sun_times = testdata.sun_position + np.testing.assert_almost_equal(sun_pos, np.array([[1000, 0, 0], [-1000, 0, 1000]])) + np.testing.assert_almost_equal(sun_vel, np.array([[-500, 500, -500], [-500, -500, 500]])) + np.testing.assert_equal(sun_times, np.array([2, 4])) + +def test_inst_position_cache(testdata): + testdata._inst_pointing_table = { + 'J2000Q0' : [1, 1], + 'J2000Q1' : [0, 0], + 'J2000Q2' : [0, 0], + 'J2000Q3' : [0, 0], + 'ET' : [0, 1], + 'TimeDependentFrames' : [-1000, -100, 1]} + testdata._body_orientation_table = { + 'J2000Q0' : [1, 0.5], + 'J2000Q1' : [0, 0.5], + 'J2000Q2' : [0, 0.5], + 'J2000Q3' : [0, 0.5], + 'ET' : [0, 1], + 'TimeDependentFrames' : [80, 1]} + testdata._inst_position_table = { + 'J2000X' : [1, 0], + 'J2000Y' : [0, 1], + 'J2000Z' : [0, 0], + 'J2000XV' : [-1, 0], + 'J2000YV' : [0, -1], + 'J2000ZV' : [0, 0], + 'ET' : [0, 1]} + sensor_pos, sensor_vel, sensor_times = testdata.sensor_position + np.testing.assert_almost_equal(sensor_pos, np.array([[1000, 0, 0], [0, 0, 1000]])) + np.testing.assert_almost_equal(sensor_vel, np.array([[-1000, 0, 0], [0, 0, -1000]])) + np.testing.assert_equal(sensor_times, np.array([0, 1])) diff --git a/tests/pytests/test_rotation.py b/tests/pytests/test_rotation.py index 2a54775..f1b0254 100644 --- a/tests/pytests/test_rotation.py +++ b/tests/pytests/test_rotation.py @@ -102,6 +102,15 @@ def test_from_euler(): assert rot.source == 0 assert rot.dest == 1 +def test_from_euler_degrees(): + rad_angles = [[np.pi/2, np.pi/2, 0], + [-np.pi/2, -np.pi/2, 0]] + degree_angles = [[90, 90, 0], + [-90, -90, 0]] + rad_rot = TimeDependentRotation.from_euler('XYZ', rad_angles, [0, 1], 0, 1) + degree_rot = TimeDependentRotation.from_euler('XYZ', degree_angles, [0, 1], 0, 1, degrees=True) + np.testing.assert_almost_equal(rad_rot.quats, degree_rot.quats) + def test_from_matrix(): mat = [[0, 0, 1], [1, 0 ,0], diff --git a/tests/pytests/test_transformation.py b/tests/pytests/test_transformation.py index 334c0f1..0b49c08 100644 --- a/tests/pytests/test_transformation.py +++ b/tests/pytests/test_transformation.py @@ -2,7 +2,54 @@ import pytest import numpy as np from ale.rotation import ConstantRotation, TimeDependentRotation -from ale.transformation import FrameChain +from ale.transformation import FrameChain, create_rotations + +def test_create_quaternion_rotations(): + test_table = { + 'J2000Q0' : [1, 0.5], + 'J2000Q1' : [0, 0.5], + 'J2000Q2' : [0, 0.5], + 'J2000Q3' : [0, 0.5], + 'ET' : [0, 1], + 'TimeDependentFrames' : [-1000, -100, 1]} + rotations = create_rotations(test_table) + assert len(rotations) == 1 + assert rotations[0].source == 1 + assert rotations[0].dest == -1000 + np.testing.assert_equal(rotations[0].times, np.array([0, 1])) + np.testing.assert_almost_equal(rotations[0].quats, np.array([[0, 0, 0, 1], [0.5, 0.5, 0.5, 0.5]])) + +def test_create_two_rotations(): + test_table = { + 'J2000Q0' : [1, 0.5], + 'J2000Q1' : [0, 0.5], + 'J2000Q2' : [0, 0.5], + 'J2000Q3' : [0, 0.5], + 'ET' : [0, 1], + 'TimeDependentFrames' : [-1000, -100, 1], + 'ConstantRotation' : [0, 0, 1, 1, 0 , 0, 0, 1, 0], + 'ConstantFrames' : [-1020, -1000]} + rotations = create_rotations(test_table) + assert len(rotations) == 2 + assert rotations[1].source == -1000 + assert rotations[1].dest == -1020 + np.testing.assert_almost_equal(rotations[1].quat, np.array([0.5, 0.5, 0.5, 0.5])) + +def test_create_euler_rotations(): + test_table = { + 'J2000Ang1' : [0, 0, 10], + 'J2000Ang2' : [90, -90, 2], + 'J2000Ang3' : [90, -90, 1], + 'CkTableStartTime' : 10, + 'CkTableEndTime' : 12, + 'CkTableOriginalSize' : 2, + 'TimeDependentFrames' : [-1000, -100, 1]} + rotations = create_rotations(test_table) + assert len(rotations) == 1 + assert rotations[0].source == 1 + assert rotations[0].dest == -1000 + np.testing.assert_equal(rotations[0].times, np.array([10, 12])) + np.testing.assert_almost_equal(rotations[0].quats, np.array([[0.5, 0.5, 0.5, 0.5], [0, 0, 0, 1]])) @pytest.fixture(scope='function') def frame_tree(request): @@ -24,11 +71,8 @@ def frame_tree(request): ConstantRotation(np.array([1.0/np.sqrt(2), 0, 0, 1.0/np.sqrt(2)]), 3, 2), ConstantRotation(np.array([1.0/np.sqrt(2), 0, 0, 1.0/np.sqrt(2)]), 4, 1) ] - frame_chain.add_edge(1, 2, rotation = rotations[0].inverse()) frame_chain.add_edge(2, 1, rotation = rotations[0]) - frame_chain.add_edge(2, 3, rotation = rotations[1].inverse()) frame_chain.add_edge(3, 2, rotation = rotations[1]) - frame_chain.add_edge(1, 4, rotation = rotations[2].inverse()) frame_chain.add_edge(4, 1, rotation = rotations[2]) return frame_chain, rotations @@ -112,13 +156,9 @@ def test_last_time_dependent_frame_between(): np.array([1]), 4, 1), ConstantRotation(np.array([1.0/np.sqrt(2), 0, 0, 1.0/np.sqrt(2)]), 5, 4) ] - frame_chain.add_edge(1, 2, rotation = rotations[0].inverse()) frame_chain.add_edge(2, 1, rotation = rotations[0]) - frame_chain.add_edge(2, 3, rotation = rotations[1].inverse()) frame_chain.add_edge(3, 2, rotation = rotations[1]) - frame_chain.add_edge(1, 4, rotation = rotations[2].inverse()) frame_chain.add_edge(4, 1, rotation = rotations[2]) - frame_chain.add_edge(4, 5, rotation = rotations[3].inverse()) frame_chain.add_edge(5, 4, rotation = rotations[3]) # last frame from node 1 to node 3 -- GitLab