updates data_isis to return positions in the correct frame (#214) (c44da0bd) · Commits · aflab / astrogeology / Ale

ale/base/data_isis.py

+141 −205

Original line number	Diff line number	Diff line
		@@ -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
		Parse an ISIS table into a dictionary.

		Parameters
		----------
		field_label : PVLModule
		The field label
		table_label : PVLModule
		The ISIS table label
		data : bytes
		The binary component of the ISIS table

		Returns
		-------
		int :
		The size of the one entry in bytes
		dict :
		The table as a dictionary with the keywords from the label and the
		binary data
		"""
		data_sizes = {
		'Integer' : 4,
		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',
		'Text' : 1}
		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
		'Real' : 'f'}

		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.

		Parameters
		----------
		table_label : PVLModule
		The table label

		data : bytes
		The binary data for the entire table

		Returns
		-------
		dict :
		The table as a dict
		"""
		# 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']] )

		return results
		: 2darray
		Array of rotated positions
		: array
		Array of times for the positions

		"""
		# 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]

ale/encoders.py

+2 −0

Original line number	Diff line number	Diff line
		@@ -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):

ale/rotation.py

+14 −2

Original line number	Diff line number	Diff line
		@@ -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)

ale/transformation.py

+79 −0

Original line number	Diff line number	Diff line
		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)

tests/pytests/test_data_isis.py

+122 −1

File changed.

Preview size limit exceeded, changes collapsed.