r''' VEX spice related module
Tutorial is available at :ref:`tutorial_vex_spice`.
Following the tutorial, users have to download the kernels and
:file:`.irfpyrc` file should have following two entries
to specify the location of the kernel.
::
[vexpvat]
spiceroot = /Volumes/scidata/data/venus/spice/kernels/
To use the vexspice module, user has to initialized via :func:`init` function.
>>> init()
Then use :func:`get_position` function to retrieve the position of Venus Express.
>>> print(get_position(datetime.datetime(2008, 12, 5, 13, 0, 0)))
[ -5700.13781279 19446.89874591 -41646.92614695]
.. note::
Since March 2017, ESA looks to provide MK file, which contains the
meta data of SPICE. The current version just ignores this file.
For the future development, the MK file will be used.
'''
import os
import logging
import datetime
import warnings as _warnings
import numpy as np
import spiceypy as spice
import irfpy.util.irfpyrc
logger = logging.getLogger(__name__)
[docs]def isdb():
""" Return whether the SPICE kernels are located in the correct location specified in .irfpyrc.
"""
rc = irfpy.util.irfpyrc.Rc()
spiceroot = rc.get('vexpvat', 'spiceroot')
if spiceroot is None:
return False
if not os.path.isdir(spiceroot):
return False
return True
[docs]class VexSpiceFactory:
def __init__(self):
self.rc = irfpy.util.irfpyrc.Rc()
self.loaded = []
self.exclude_kernel = []
[docs] def showstatus(self):
print("Loaded kernels = %d" % len(self.loaded))
for k in self.loaded:
print(' -- %s' % k)
[docs] def get_filename_list(self, level='full', dataonly=True):
"""
:param level:
:param dataonly:
:return:
"""
return self.get_full(dataonly=dataonly)
[docs] def get_full(self, dataonly=True):
""" Get the kernel file name list
:param dataonly:
:return: The file name list of the SPICE kernels
"""
spiceroot = self.rc.get('vexpvat', 'spiceroot')
if spiceroot is None:
logger.warning('Specify the root of spice in .irfpyrc'
'by [vexpvat] spiceroot entry')
raise RuntimeError('')
logger.debug('spiceroot = %s' % spiceroot)
if not os.path.isdir(spiceroot):
_warnings.warn('It looks the folder specified by .irfpyrc does not exist.\n' +
'Specified folder: {}\n'.format(spiceroot) +
'Check the file ~/.irfpyrc or ~/irfpy.rc')
logging.debug('spiceroot= %s' % spiceroot)
kfiles = []
for p, d, f in os.walk(spiceroot):
for ff in f:
# print ff
if ff.endswith('.txt'):
continue
if ff.lower().endswith('.tm'):
continue
if ff.endswith('~'):
continue
logging.debug('Adding %s into kernel list' % ff)
kfiles.append(os.path.join(p, ff))
for fname in self.exclude_kernel:
logger.debug('Remove {}'.format(fname))
kfiles.remove(fname)
return kfiles
[docs] def furnsheach(self, fname):
fabsname = os.path.abspath(fname)
if not fabsname in self.loaded:
if not os.path.exists(fabsname):
raise IOError('Kernel file {} not found'.format(fabsname))
logger.info('Furnsh: {}'.format(fabsname))
spice.furnsh(fabsname)
self.loaded.append(fabsname)
[docs] def furnsh(self):
''' Furnsh the kernels.
'''
fnames = self.get_full()
for fn in fnames:
try:
self.furnsheach(fn)
except spice.stypes.SpiceyError:
logger.warning("Failed to load {}. Usually no problem. Let's just continue.".format(fn))
_factory = VexSpiceFactory()
[docs]def isinitialized():
_warnings.warn('{}.isinitialized() deprecated. Do nothing.'.format(__name__))
[docs]def init(update=True):
''' Initialize the VEX spice kernels
:keyword update: Not used anymore.
'''
_factory.furnsh()
[docs]def init_noatt(update=True):
''' Furnsh default kernels from the irfpyrc's mk file
:param update: Update the mk file, but not for attitude file.
.. warning::
This function is now the same as :func:`init`.
'''
init(update=update)
[docs]def furnsh(kernel_file):
''' Furnsh additional kernel file
'''
spice.furnsh(kernel_file)
[docs]def updatemk(noatt=False, mkfile=None):
''' Update the mk file.
Walk throught the path under ``spiceroot``,
create ``mk`` file.
Excluding ``*.txt``
>>> updatemk()
'''
_warnings.warn("{}.updatemk() is deprecated. Do nothing.".format(__name__))
[docs]def get_position(t, target='VEX', origin='VENUS', frame='VSO', correction='LT+S'):
''' Get the position. By default, VEX in VSO.
>>> init()
>>> pos = get_position(datetime.datetime(2006, 12, 2, 15, 30))
>>> print('%.2f %.2f %.2f' % (pos[0], pos[1], pos[2])) # may change slightly
-985.40 -2387.85 -68309.44
'''
return get_posvel(t, target=target, origin=origin, frame=frame, correction=correction)[:3]
[docs]def get_velocity(t, target='VEX', origin='VENUS', frame='VSO', correction='LT+S'):
''' Get the velocity. By default, VEX in VSO.
>>> init()
>>> vel = get_velocity(datetime.datetime(2006, 12, 2, 15, 30))
>>> print('%.2f %.2f %.2f' % (vel[0], vel[1], vel[2])) # may change slightly
0.85 0.26 -0.70
'''
return get_posvel(t, target=target, origin=origin, frame=frame, correction=correction)[3:]
[docs]def get_positions(tlist, target='VEX', origin='Venus', frame='VSO', correction='LT+S'):
''' Get the time series of the position.
:param tlist: Iterable of ``datetime.datetime`` object.
:returns: (N, 3) array of the position vector.
>>> init()
>>> tlist = [datetime.datetime(2012, 1, 1, 12), datetime.datetime(2012, 1, 1, 13)]
>>> positions = get_positions(tlist)
>>> print(positions.shape)
(2, 3)
>>> print("{p[0]:.2f} {p[1]:.2f} {p[2]:.2f}".format(p=positions[0, :]))
-7673.50 26101.15 -37544.75
>>> print("{p[0]:.2f} {p[1]:.2f} {p[2]:.2f}".format(p=positions[1, :]))
-7858.35 25955.30 -45445.77
'''
poss = [get_position(t, target=target, origin=origin, frame=frame, correction=correction) for t in tlist]
return np.array(poss)
[docs]def get_velocities(tlist, target='VEX', origin='Venus', frame='VSO', correction='LT+S'):
''' Get the time series of the velocity.
:param tlist: Iterable of ``datetime.datetime`` object.
:returns: (N, 3) array of the position vector.
>>> init()
>>> tlist = [datetime.datetime(2012, 1, 1, 12), datetime.datetime(2012, 1, 1, 13)]
>>> velocities = get_velocities(tlist)
>>> print(velocities.shape)
(2, 3)
>>> print("{p[0]:.2f} {p[1]:.2f} {p[2]:.2f}".format(p=velocities[0, :]))
-0.09 0.09 -2.40
>>> print("{p[0]:.2f} {p[1]:.2f} {p[2]:.2f}".format(p=velocities[1, :]))
-0.02 -0.16 -2.00
'''
poss = [get_velocity(t, target=target, origin=origin, frame=frame, correction=correction) for t in tlist]
return np.array(poss)
[docs]def get_posvel(t, target='VEX', origin='VENUS', correction='LT+S', frame='VSO'):
et0 = spice.str2et(t.strftime('%FT%T.%f'))
posvel, lt = spice.spkezr(target, et0, frame, correction, origin)
return np.array(posvel)
[docs]def get_scannerangle(t):
''' Return the scanner angle in degrees.
'''
### The scanner angle is the angle between +Xsc and +Xsaf.
et0 = spice.str2et(t.strftime('%FT%T.%f'))
sc2saf = spice.pxform('VEX_SPACECRAFT', 'VEX_ASPERA4_SAF', et0)
saf = spice.mxv(sc2saf, [1, 0, 0])
sdeg = np.rad2deg(np.arctan2(saf[1], saf[0]))
return sdeg
[docs]def get_lonlatalt(t, target='VEX', body='VENUS', frame='IAU_VENUS', re=6052, f=0):
"""
codeauthor:: Moa Persson
:param t: input specific time in datetime format
:param target: target body of which longitude we seek, standard is Venus Express
:param body: body of which longitude we look at
:param frame: has to be body fixed, standard is IAU Venus
:param re: radius of body, standard Venus radius 6052 km
:param f: Flattening coefficient, :math:`(re-rp) / re`
where :math:`rp` is the polar radius of the spheroid, and the
units of :math:`rp` match those of :math:`re`. Standard for Venus: zero oblateness.
:returns: longitude, latitude and altitude of target. longitude and latitude as radians.
altitude as the point above reference spheloid.
.. seealso:: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/recpgr_c.html
>>> init()
>>> t = datetime.datetime(2006, 12, 2, 15, 30)
>>> lon, lat, alt = get_lonlatalt(t)
>>> print('%.2f %.2f %.2f' % (lon, lat, alt))
1.56 -1.56 62306.26
"""
pos_x, pos_y, pos_z = get_position(t, frame=frame, origin=body, target=target)
lon, lat, alt = spice.recpgr(body=body, rectan=[pos_x, pos_y, pos_z], re=re, f=f)
return lon, lat, alt
[docs]def get_lonlat(*args, **kwds):
""" Alias to :func:`get_lonlatalt`.
"""
return get_lonlatalt(*args, **kwds)
[docs]def get_lst(t, lon, bodyID=299):
"""
codeauthor:: Moa Persson
:param t: input specific time in datetime format
:param lon: input longitude on the surface
:param bodyID: spice ID of body, standard 299 is Venus
:returns: local solar time float and string
>>> init()
>>> t = datetime.datetime(2006, 12, 2, 15, 30)
>>> lon, lat, alt = get_lonlat(t)
>>> lst_frac, time = get_lst(t, lon)
>>> print('%.2f' % (lst_frac), time) # may change slightly
10.84 10:50:17
"""
et0 = spice.str2et(t.strftime('%FT%T.%f'))
hr, mn, sc, time, ampm = spice.et2lst(et0, body=bodyID, lon=lon, typein='PLANETOGRAPHIC', timlen=256, ampmlen=256)
lst_frac = hr + mn / 60 + sc / 3600
return lst_frac, time
[docs]def get_sza(t, origin='Venus', target='VEX'):
"""
codeauthor:: Moa Persson
:param t: input specific time in datetime format
:param origin: origin of the coordinate system vector
:param target: target of the coordinate system vector
:returns: solar zenith angle (sza) of the target in the frame of the origin body at the specified t
>>> init()
>>> sza = get_sza(datetime.datetime(2006, 12, 2, 15, 30))
>>> print('%.2f' % (sza))
90.83
"""
# Position of target and Sun in frame
pos_x, pos_y, pos_z = get_position(t, origin=origin, target=target)
pos_x_sun, pos_y_sun, pos_z_sun = get_position(t, origin=origin, target='Sun')
# Length of vectors
altitude = np.sqrt(pos_x ** 2 + pos_y ** 2 + pos_z ** 2) # Altitude above the reference spheroid
alt_sun = np.sqrt(pos_x_sun ** 2 + pos_y_sun ** 2 + pos_z_sun ** 2) # Distance to the sun
# Normalise vectors
pos_vex = np.array([pos_x, pos_y, pos_z]) / altitude
pos_sun = np.array([pos_x_sun, pos_y_sun, pos_z_sun]) / alt_sun
# Calculate Sun-Venus-Spacecraft angle, i.e. SZA
svsc = 180 / np.pi * np.arccos(np.dot(pos_sun, pos_vex))
return svsc
[docs]def convert_matrix(t, fromframe="VEX_SPACECRAFT", toframe="VSO"):
''' Return conversion matrix as a numpy array.
:param t: ``datatime.datetime`` object specifying time.
:returns: Numpy array of the conversion matrix. ``numpy.pxform`` is used.
Conversion between IMA frame and SC frame is
>>> from irfpy.vexpvat import vexspice
>>> vexspice.init()
>>> import datetime
>>> t = datetime.datetime.now() # This is just a placeholder, as they are fixed each other.
>>> matx = convert_matrix(t, fromframe='VEX_ASPERA4_IMA', toframe='VEX_SPACECRAFT')
>>> from irfpy.util.with_context import printoptions
>>> with printoptions(suppress=True, precision=3): print(matx)
[[-0. 1. -0.]
[ 1. 0. 0.]
[ 0. -0. -1.]]
'''
# frame0 = spicename(fromframe)
# frame1 = spicename(toframe)
et0 = spice.str2et(t.strftime('%FT%T.%f'))
return np.array(spice.pxform(fromframe, toframe, et0))