In [332]:
import json
import math
import numpy as np
import pvl
import spiceypy as spice
In [333]:
# 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 [334]:
# Furnish the IK Kernel
# spice.furnsh("../../tests/data/msgr_mdis_v160.ti")
# ikid = 236820
# spice.furnsh('../../tests/data/msgr_v231.tf') # Frames Kernel mapping between frames
# spice.furnsh('../../tests/data/pck00010_msgr_v23.tpc') # PC Kernel with planetary attitude and body information
# spice.furnsh('../../tests/data/msgr_dyn_v600.tf')
# spice.furnsh('../../tests/data/msgr_de405_de423s.bsp')
# spice.furnsh('../../tests/data/msgr_040803_150430_150430_od431sc_2.bsp')
# spice.furnsh('../../tests/data/msgr_mdis_sc050727_100302_sub_v1.bc')
# spice.furnsh('../../tests/data/msgr_mdis_gm040819_150430v1.bc') # MDIS Instrument Pointing (236890)
# spice.furnsh('../../tests/data/msgr_1304_v02.bc') #Bus and Spacecraft Panels position
# spice.furnsh('../../tests/data/naif0011.tls')
# spice.furnsh('../../tests/data/messenger_2548.tsc')
ikid = 236800
msgr = "/usgs/cpkgs/isis3/data/messenger/"
# Load kernels same order ISIS Spice::init() does
# Frame
# TargetAttitudeShape
spice.furnsh(msgr+"kernels/pck/pck00010_msgr_v23.tpc")
# Instrument
spice.furnsh(msgr+"kernels/ik/msgr_mdis_v160.ti")
# InstrumentAddendum
spice.furnsh(msgr+"kernels/iak/mdisAddendum009.ti")
# LeapSecond
spice.furnsh("/usgs/cpkgs/isis3/data/base/kernels/lsk/naif0012.tls")
# SpacecraftClock
spice.furnsh(msgr+"kernels/sclk/messenger_2548.tsc")
# Extra
# TargetPosition
spice.furnsh(msgr+"kernels/tspk/de423s.bsp")
# InstrumentPointing
spice.furnsh(msgr+"kernels/ck/msgr20150409.bc")
spice.furnsh(msgr+"kernels/ck/msgr20150410.bc")
spice.furnsh(msgr+"kernels/ck/msgr20150411.bc")
spice.furnsh(msgr+"kernels/ck/msgr20150412.bc")
spice.furnsh(msgr+"kernels/ck/msgr20150413.bc")
spice.furnsh(msgr+"kernels/ck/msgr20150414.bc")
spice.furnsh(msgr+"kernels/ck/msgr20150415.bc")
spice.furnsh(msgr+"kernels/ck/msgr20150416.bc")
spice.furnsh(msgr+"kernels/ck/1072683119_1965_mdis_atthist.bc")
spice.furnsh(msgr+"kernels/ck/1072716050_291010_mdis_pivot_pvtres.bc")
spice.furnsh(msgr+"kernels/fk/msgr_v231.tf")
# InstrumentPosition
spice.furnsh(msgr+"kernels/spk/msgr_20040803_20150430_od431sc_2.bsp")
In [335]:
# 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)
What is the minimum information needed to convert to and from the MDIS-NAC frame?
236: The code for the messenger spacecraft 4: Indicates that the frame is a TK Frame 236820: The code for the MDIS-NAC Frame
In [336]:
# Type 4 is TK - the mdis-nac is constant w.r.t. another frame
print(spice.frinfo(-236820))
# Verify that this is messenger
print(spice.bods2c('MESSENGER'))
In [337]:
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 [338]:
# Load the ISIS Cube header
# header = pvl.load('../../tests/data/EN1007907102M.cub')
header = pvl.load("../../tests/data/CW1071364100B_IU_5.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')
filter_num = find_in_dict(header, 'OriginalFilterNumber')
ikid += filter_num
# 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 [339]:
# 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
This is almost identical to ISIS3. Verified using campt.
In [340]:
# This needs to be sensor origin - like the comment below says...
# This is almost there - position w.r.t. the spacecraft, but not the camera yet.
loc_j2000, ltc = spice.spkpos(isd['target_name'], isd['ephemeris_time'], 'J2000', 'LT+S', 'MESSENGER')
loc_j2000
print("Spacecraft in J2000:\n\t{0:.15f}\n\t{1:.15f}\n\t{2:.15f}".format(loc_j2000[0], loc_j2000[1], loc_j2000[2]))
j2000_to_mercury = spice.pxform('J2000', 'IAU_MERCURY', isd['ephemeris_time'])
loc_mercury = spice.mxv(j2000_to_mercury, loc_j2000)
print("Spacecraft in MERCURY:\n\t{0:.15f}\n\t{1:.15f}\n\t{2:.15f}".format(loc_mercury[0], loc_mercury[1], loc_mercury[2]))
# Direct to Mercury body-fixed position is slightly (neglibly) off (campt SpacecraftPosition)
loc_direct, lt_direct = spice.spkpos(isd['target_name'], isd['ephemeris_time'], 'IAU_MERCURY', 'LT+S', 'MESSENGER')
print("Spacecraft direct spkpos in Mercury frame:")
loc_direct
print("\t{0:.15f}\n\t{1:.15f}\n\t{2:.15f}".format(loc_direct[0], loc_direct[1], loc_direct[2]))
# Try to compare with the spkez routine
#loc_ez, ltc_ez = spice.spkez(199, isd['ephemeris_time'], 'IAU_MERCURY', 'LT+S', -236)
#print("{0:.15f} \n{1:.15f} \n{2:.15f}".format(loc_ez[0], loc_ez[1], loc_ez[2]))
In [351]:
loc, ltc = spice.spkpos(isd['target_name'], isd['ephemeris_time'], 'J2000', 'LT+S', 'MESSENGER')
# Transform from J2000 to IAU_MERCURY (body-fixed) frame
rotation = spice.pxform('J2000', 'IAU_MERCURY', isd['ephemeris_time'])
loc = spice.mxv(rotation, loc)
isd['x_sensor_origin'] = loc[0] * -1000
isd['y_sensor_origin'] = loc[1] * -1000
isd['z_sensor_origin'] = loc[2] * -1000
In [342]:
print("{0:.15f}".format(isd['x_sensor_origin']))
print("{0:.15f}".format(isd['y_sensor_origin']))
print("{0:.15f}".format(isd['z_sensor_origin']))
In [343]:
# Get the camera position in lat, lon, height
# lat/lon are in raddians
_, radii = spice.bodvrd("MERCURY", "RADII", 3)
radii *= 1000
flattening = (radii[0] - radii[2]) / radii[0]
lon, lat, height = spice.recgeo(loc * 1000, radii[0], flattening)
# Can I use X, Y, Z body fixed in meters, or do I need to use lat, lon, height
# in rectangular coordinates?
#isd['x_sensor_origin'] = lon
#isd['y_sensor_origin'] = lat
#isd['z_sensor_origin'] = height
In [344]:
# This needs to go from object space to image space - right now, it goes from object space (Mercury)
# to sensor position space (MDIS_NAC)
# camera2bodyfixed = spice.pxform('MSGR_MDIS_NAC', 'IAU_MERCURY', isd['ephemeris_time'])
camera2bodyfixed = spice.pxform('MSGR_MDIS_WAC', 'IAU_MERCURY', isd['ephemeris_time'])
camopk = spice.m2eul(camera2bodyfixed, 3, 2, 1)
isd['omega'] = camopk[2]
isd['phi'] = camopk[1]
isd['kappa'] = camopk[0]
print(camopk)
# Rotation around Y, Z from camera to focal plane
# focal = np.zeros((3,3))
# focal[0][0] = 1
# focal[1][1] = -1
# focal[2][2] = -1
# focalrot = spice.mxmt(focal, camera2bodyfixed) # np.dot().T
# # Get the length of the vector
# xyzlength = loc[0]**2 + loc[1]**2
# xylength = math.sqrt(xyzlength)
# xyzlength = math.sqrt(xyzlength + loc[2]**2)
# slon = loc[1] / xylength
# clon = loc[0] / xylength
# slat = loc[2] / xyzlength
# clat = xylength / xyzlength
# #Ocentric to ographic rotation
# oo = np.zeros((3,3))
# oo[:,0] = -slon, clon, 0.0
# oo[:,1] = -slat * clon, -slat * slon, clat
# oo[:,2] = clat*clon, clat*slon, slat
# rot = spice.mxm(focalrot, oo) # np.dot
# print(rot)
# opk = spice.m2eul(rot, 3, 2, 1)
# print(opk)
# print(list(map(math.degrees, opk)))
#isd['omega'] = opk[2]
#isd['phi'] = opk[1]
#isd['kappa'] = opk[0]
In [345]:
camera2bodyfixed
Out[345]:
In [ ]:
# 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 [346]:
#print(list(map(degrees,isd['omega'], isd['phi'], isd['kappa'])))
In [347]:
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 [348]:
!cat isd.isd
In [349]:
# 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 [350]:
# OPK
isd['x_sensor_origin'] =
isd['y_sensor_origin'] =
isd['z_sensor_origin'] =
isd['omega'] =
isd['phi'] =
isd['kappa'] =
In [ ]:
# ISD Search Information - totally fabricated.
isd['min_elevation'] = -1.0
isd['max_elevation'] = 1.0
In [ ]:
In [ ]:
isd
In [ ]: