In [55]:
import json
import math
import numpy as np
import pvl
import spiceypy as spice
In [56]:
# Utility Func
def find_in_dict(obj, key):
"""
Recursively find an entry in a dictionary
Parameters
----------
obj : dict
The dictionary to search
key : str
The key to find in the dictionary
Returns
-------
item : obj
The value from the dictionary
"""
if key in obj:
return obj[key]
for k, v in obj.items():
if isinstance(v,dict):
item = find_in_dict(v, key)
if item is not None:
return item
In [57]:
ikid = 236820
# Load kernels same order ISIS Spice::init() does
# Frame
# TargetAttitudeShape
spice.furnsh('../../tests/data/pck00010_msgr_v23.tpc')
# Instrument
spice.furnsh("../../tests/data/msgr_mdis_v160.ti")
# InstrumentAddendum
spice.furnsh('../../tests/data/mdisAddendum009.ti')
# LeapSecond
spice.furnsh('../../tests/data/naif0012.tls')
# SpacecraftClock
spice.furnsh('../../tests/data/messenger_2548.tsc')
# Extra
# TargetPosition
spice.furnsh('../../tests/data/de423s.bsp')
# InstrumentPointing
spice.furnsh('../../tests/data/msgr20130404.bc')
spice.furnsh('../../tests/data/msgr20130405.bc')
spice.furnsh('../../tests/data/msgr20130406.bc')
spice.furnsh('../../tests/data/msgr20130407.bc')
spice.furnsh('../../tests/data/msgr20130408.bc')
spice.furnsh('../../tests/data/msgr20130409.bc')
spice.furnsh('../../tests/data/msgr20130410.bc')
spice.furnsh('../../tests/data/msgr20130411.bc')
spice.furnsh('../../tests/data/1072683119_1965_mdis_atthist.bc')
spice.furnsh('../../tests/data/1072716050_291010_mdis_pivot_pvtres.bc')
spice.furnsh('../../tests/data/msgr_v231.tf')
# InstrumentPosition
spice.furnsh('../../tests/data/msgr_20040803_20150430_od431sc_2.bsp')
In [58]:
# Create the ISD object
isd = {}
# Load information from the IK kernel
isd['focal_length'] = spice.gdpool('INS-{}_FOCAL_LENGTH'.format(ikid), 0, 1)
isd['focal_length_epsilon'] = spice.gdpool('INS-{}_FL_UNCERTAINTY'.format(ikid), 0, 1)
isd['nlines'] = spice.gipool('INS-{}_PIXEL_LINES'.format(ikid), 0, 1)
isd['nsamples'] = spice.gipool('INS-{}_PIXEL_SAMPLES'.format(ikid), 0, 1)
isd['original_half_lines'] = isd['nlines'] / 2.0
isd['original_half_samples'] = isd['nsamples'] / 2.0
isd['pixel_pitch'] = spice.gdpool('INS-{}_PIXEL_PITCH'.format(ikid), 0, 1)
isd['ccd_center'] = spice.gdpool('INS-{}_CCD_CENTER'.format(ikid), 0, 2)
isd['ifov'] = spice.gdpool('INS-{}_IFOV'.format(ikid), 0, 1)
isd['boresight'] = spice.gdpool('INS-{}_BORESIGHT'.format(ikid), 0, 3)
isd['transx'] = spice.gdpool('INS-{}_TRANSX'.format(ikid), 0, 3)
isd['transy'] = spice.gdpool('INS-{}_TRANSY'.format(ikid), 0, 3)
isd['itrans_sample'] = spice.gdpool('INS-{}_ITRANSS'.format(ikid), 0, 3)
isd['itrans_line'] = spice.gdpool('INS-{}_ITRANSL'.format(ikid), 0, 3)
isd['odt_x'] = spice.gdpool('INS-{}_OD_T_X'.format(ikid), 0, 10)
isd['odt_y'] = spice.gdpool('INS-{}_OD_T_Y'.format(ikid), 0, 10)
isd['starting_detector_sample'] = spice.gdpool('INS-{}_FPUBIN_START_SAMPLE'.format(ikid), 0, 1)
isd['starting_detector_line'] = spice.gdpool('INS-{}_FPUBIN_START_LINE'.format(ikid), 0, 1)
In [59]:
def distort_focal_length(coeffs, t):
"""
Compute the distorted focal length
Parameters
----------
coeffs : iterable
of coefficient values
t : float
temperature in C
Returns
-------
focal_length : float
the temperature adjusted focal length
"""
focal_length = coeffs[0]
for i in range(1, len(coeffs[1:])):
focal_length += coeffs[i]*t**i
return focal_length
In [60]:
# Load the ISIS Cube header
header = pvl.load('../../tests/data/EN1007907102M.cub')
isd['instrument_id'] = find_in_dict(header, 'InstrumentId')
isd['spacecraft_name'] = find_in_dict(header, 'SpacecraftName')
isd['target_name'] = find_in_dict(header, 'TargetName')
# Get the radii from SPICE
rad = spice.bodvrd(isd['target_name'], 'RADII', 3)
radii = rad[1]
isd['semi_major_axis'] = rad[1][0]
isd['semi_minor_axis'] = rad[1][1]
# Get temperature from SPICE and adjust focal length
spice.gdpool('INS-{}_FOCAL_LENGTH'.format(ikid), 0, 1)
temp_coeffs = spice.gdpool('INS-{}_FL_TEMP_COEFFS'.format(ikid), 0, 6)
temp = find_in_dict(header, 'FocalPlaneTemperature').value
isd['focal_length'] = distort_focal_length(temp_coeffs, temp)
Now we need to compute time. - This has been verified as correct using campt.
In [61]:
# Here convert the sclock
sclock = find_in_dict(header, 'SpacecraftClockCount')
exposure_duration = find_in_dict(header, 'ExposureDuration')
exposure_duration = exposure_duration.value * 0.001 # Scale to seconds
# Get the instrument id, and, since this is a framer, set the time to the middle of the exposure
spacecraft_id = spice.bods2c('MESSENGER')
et = spice.scs2e(spacecraft_id, sclock)
et += (exposure_duration / 2.0)
isd['ephemeris_time'] = et
In [62]:
# Comparing spacecraft position direct to Mercury body-fixed frame with J2000 to Mercury body-fixed
loc_direct, _ = spice.spkpos(isd['target_name'], isd['ephemeris_time'], 'IAU_MERCURY', 'LT+S', 'MESSENGER')
print(loc_direct * -1000)
# Target=Mercury, Observer=Messenger, go to J2000 first (to match ISIS3)
loc, ltc = spice.spkpos(isd['target_name'], isd['ephemeris_time'], 'J2000', 'LT+S', 'MESSENGER')
# Get the transformation from J2000 to Mercury body-fixed frame
rotation = spice.pxform('J2000', 'IAU_MERCURY', isd['ephemeris_time'])
loc = spice.mxv(rotation, loc)
# Scale to meters and reverse direction (since target=Mercury and observer=messenger)
# we want vector from Mercury origin to spacecraft position
isd['x_sensor_origin'] = loc[0] * -1000
isd['y_sensor_origin'] = loc[1] * -1000
isd['z_sensor_origin'] = loc[2] * -1000
In [63]:
print(isd['x_sensor_origin'])
print(isd['y_sensor_origin'])
print(isd['z_sensor_origin'])
In [64]:
# Get the rotation angles from MDIS NAC frame to Mercury body-fixed frame
camera2bodyfixed = spice.pxform('MSGR_MDIS_NAC','IAU_MERCURY', isd['ephemeris_time'])
opk = spice.m2eul(camera2bodyfixed, 3, 2, 1)
isd['omega'] = opk[2]
isd['phi'] = opk[1]
isd['kappa'] = opk[0]
In [65]:
camera2bodyfixed
Out[65]:
In [70]:
# Get the sun's position relative to a Mercury-fixed frame.
target = "SUN"
et = isd['ephemeris_time']
reference_frame = "IAU_MERCURY"
light_time_correction = "LT+S"
observer = "MERCURY"
sun_state, lt = spice.spkezr(target,
et,
reference_frame,
light_time_correction,
observer)
# Convert to meters
isd['x_sun_position'] = sun_state[0] * 1000
isd['y_sun_position'] = sun_state[1] * 1000
isd['z_sun_position'] = sun_state[2] * 1000
print("SUN POSITION (m): {} {} {}".format(sun_state[0]*1000,
sun_state[1]*1000,
sun_state[2]*1000))
# Get velocity of mdis nac (right now it is getting velocity of spacecraft)
target = "MESSENGER"
et = isd['ephemeris_time']
reference_frame = "IAU_MERCURY"
light_time_correction = "LT+S"
observer = "MERCURY"
messenger_state, lt = spice.spkezr(target,
et,
reference_frame,
light_time_correction,
observer)
print("MESSENGER VELOCITY (m/s): {} {} {}".format(messenger_state[3]*1000,
messenger_state[4]*1000,
messenger_state[5]*1000))
v,_ = spice.spkezr(observer, et, reference_frame, light_time_correction, target)
print(v)
In [67]:
class NumpyAwareJSONEncoder(json.JSONEncoder):
def default(self, obj):
if isinstance(obj, np.ndarray) and obj.ndim == 1:
lobj = obj.tolist()
if len(lobj) == 1:
return lobj[0]
else:
return lobj
return json.JSONEncoder.default(self, obj)
with open('isd.isd', 'w') as f:
f.write(json.dumps(isd, f, cls=NumpyAwareJSONEncoder, sort_keys=True, indent=4))
In [68]:
!cat isd.isd
In [ ]:
In [69]:
%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
from itertools import product, combinations
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_aspect('equal')
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = 1 * np.outer(np.cos(u), np.sin(v))
y = 1 * np.outer(np.sin(u), np.sin(v))
z = 1 * np.outer(np.ones(np.size(u)), np.cos(v))
#for i in range(2):
# ax.plot_surface(x+random.randint(-5,5), y+random.randint(-5,5), z+random.randint(-5,5), rstride=4, cstride=4, color='b', linewidth=0, alpha=0.5)
elev = 10.0
rot = 80.0 / 180 * np.pi
ax.plot_surface(x, y, z, rstride=4, cstride=4, color='b', linewidth=0, alpha=0.5)
# Plot the origin of the body fixed coordinate system
ax.scatter(0, 0, 0, 'ko', s=50)
# Plot the spacecraft position, normalized to the radius of the sphere
x = isd['x_sensor_origin'] / radii[0]
y = isd['y_sensor_origin'] / radii[0]
z = isd['z_sensor_origin'] / radii[0]
delta = 0.075
px = float(x + delta)
py = float(y + delta)
pz = float(z + delta)
nx = float(x - delta)
ny = float(y - delta)
nz = float(z - delta)
coords = np.array([(px, py, pz),
(px, py, nz),
(px, ny, nz),
(px, ny, pz),
(px, py, pz),
(nx, py, pz),
(nx, py, nz),
(nx, ny, nz),
(nx, ny, pz),
(nx, py, pz),
(nx, ny, pz),
(px, ny, pz),
(px, ny, nz),
(nx, ny, nz),
(nx, py, nz),
(px, py, nz)])
ax.plot_wireframe(coords[:,0], coords[:,1], coords[:,2],color= 'r')
#print(isd['boresight'])
ax.quiver(x, y, z, 0,0,1, color='green', pivot='tail') # Boresight
print(camera2bodyfixed.dot(np.array([0, 0, isd['focal_length']])))
ax.quiver(x, y, z, *camera2bodyfixed.dot(np.array([0, 0, isd['focal_length']])), color='black', pivot='tail')
In [ ]:
In [ ]:
spice.spkez(-236820, isd['ephemeris_time'], 'IAU_MERCURY', "NONE", 199 )
In [ ]:
In [ ]:
!cat isd.isd
In [14]:
# This is hard coded in Anne's script
isisFocalPlan2SocetPlate = np.eye(3)
isisFocalPlan2SocetPlate[1,1] = -1.0
isisFocalPlan2SocetPlate[2,2] = -1.0
# Grab the body fixed coordinates from SPICE
# The mercury Naif ID code is 199
nid = 199
In [15]:
# OPK
isd['x_sensor_origin'] =
isd['y_sensor_origin'] =
isd['z_sensor_origin'] =
isd['omega'] =
isd['phi'] =
isd['kappa'] =
In [16]:
# ISD Search Information - totally fabricated.
isd['min_elevation'] = -1.0
isd['max_elevation'] = 1.0
In [ ]:
In [17]:
isd
Out[17]:
In [ ]:
In [ ]: