Introduction¶
In order to follow this tutorial you need:
- spiceql (1.2.3 was used here)
- numpy
- maplotlib
Goals for this tutorial are
- Learn the flow of using SpiceQL
- Perform simple UTC to ET conversions
- Learn of the universal params
- Calculate an array of target states
- Plot those states with matplotlib
Warning
In order to use local calls (when useWeb=False
), you will need to have a local install of kernels and SpiceQL metadata generated. See the SpiceQL manual for a quick overview. We recommend anyone not familiar with SPICE kernels to use the web feature.
This tutorial is an adaptation of spiceypy's Cassini example
1. Import pyspiceql¶
2 Get start and stop times¶
We define these as utc times, but they have to be converted into ephemeris times (ETs).
useWeb=False # Set to true if you do not have local kernels
# We need our times in Ets
utc = ["Jun 20, 2004", "Dec 1, 2005"]
# this still needs LSK kernels, so set useWeb
# to True if you do not have them
et1,_ = spql.utcToEt(utc[0], useWeb=useWeb)
et2,_ = spql.utcToEt(utc[1], useWeb=useWeb)
# Create 4000 times between the start and stop
step = 4000
times = [x * (et2 - et1) / step + et1 for x in range(step)]
# print all 4000 times
print(f"Getting positions for {len(times)} ephems between {et1} and {et2}")
3. Call getTargetStates¶
pyspiceql.getTargetStates
is similar to the NAIF's SPICE toolkit's spkpos
. The main difference is in order to use spkpos
you must have kernels funished already. SpiceQL will complete the query for you under the hood to load the appropriate kernels and complete the operation. It also allows for getting results without needing local kernels.
For this reason, pyspiceql.getTargetStates
(and most other SpiceQL API calls) have additional parameters to control this behavior.
These universal params are:
- ckQualities: array of CK qualities wanted, defaults are Reconstructed and Smithed. Check the glossary for more info.
- spkQualities: array of SPK qualities wanted, defaults are Reconstructed and Smithed. Check the glossary for more info.
- useWeb: Boolean on whether or not to use the web API, set to true if you do not have local kernels
- searchKernels: Boolean on whether or not to search and furnish kernels, set to False if you are furnishing your own kernels
- fullKernelPath: Boolean, set to True to return kernels with full paths. Otherwise, you get relative paths
- limitCk: int to limit number of CKs returned, large time queries can result in many CKs
- limitSpk: int to limit the number of SPKs returned; large time queries can result in many SPKs
- kernelList: array of additional kernel paths to load on top of any returned from a query.
On limitCK/limitSPK
The defaults for these are set up for queries that are expected to span a single image. This limits SPKs to 1 which forces only the highest priority SPK to be loaded. If you are searching large ranges, set this to -1 so that you get all matching SPKs.
# All SpiceQL SPICE API calls return a pair of (data, kernels)
spqlpos, kernels = spql.getTargetStates(times,
"Cassini", # the target we wants states for
"SATURN BARYCENTER", # The observing body
"J2000", # The reference frame we want positions to
"None", # light time correction
"cassini", # SpiceQL category, this is either a mission or body,
# e.g. LRO_LROCNACL would just be LRO
useWeb=useWeb, # set to True to not use local kernels
limitSpk=-1, # -1 to get all SPKs in the time range
limitCk=-1) # -1 to get all CKs in the time range
print(f"Got {len(spqlpos)} states back")
4. Optional: Verify kernels¶
Kernels are returned as a python dict (a JSON object in C++), and come in the format kernel_type : [kernels]
. Additionally, if CK or SPK kernels were used, you will see keys in the format <body_name>_ck_quality
and <body_name>_spk_quality
. This has the highest quality acquired for SPKs and CKs respectively.
If these kernels are located locally, you can furnish them with pyspiceql.KernelSet
to integrate them further with SpiceQL or spiceypy.
5. Use matplotlib to plot results¶
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
positions = (
np.asarray(spqlpos).T
) # positions is shaped (4000, 3), let's transpose to (3, 4000) for easier indexing
fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(111, projection="3d")
ax.plot(positions[0], positions[1], positions[2])
plt.title("Cassini Position Example from Jun 20, 2004 to Dec 1, 2005")
plt.show()