# SPICE defining a topocentric frame on the Moon

I would like to define a topocentric frame for a particular spot on the moon. From the frames required reading i have this:

   \begindata

FRAME_MYTOPO             = 1234567
FRAME_1234567_NAME       = 'MYTOPO'
FRAME_1234567_CLASS      = 4
FRAME_1234567_CLASS_ID   = 1234567
FRAME_1234567_CENTER     = 399100

TKFRAME_1234567_SPEC     = 'ANGLES'
TKFRAME_1234567_RELATIVE = 'MOON_ME'
TKFRAME_1234567_ANGLES   = ( <-long>, <-colat>, 180 )
TKFRAME_1234567_AXES     = (       3,        2,   3 )
TKFRAME_1234567_UNITS    = 'DEGREES'

\begintext


However, there are two issues. What and how do I define the CENTER of the frame. That needs to be the point on the surface. I guess in the above case, it has some id 399100. Plus this frame needs to be offset from the center of MOON_ME. Does this happen automatically because of the specified center point? And how do I specify the center point? (I think basically in spherical coordinates with respect to the MOON_ME it should be given by (r, colat, long)?

I believe what you'll want to do is create an SPK kernel for the ephemeris of the point on the surface of the moon and then whatever ID you give that point should be the center of the topocentric frame. That is my understanding from this info:

"This method of definition is particularly well suited for defining topocentric frames on the surface of the Earth. For example, suppose you have an SPK (ephemeris) file that specifies the location of some surface point on the Earth, and that the SPK ID code of this point is 399100. Moreover suppose you have the geodetic co-latitude (COLAT) and longitude (LONG) measured in degrees for this point. (Note that the co-latitude is the complement of latitude: latitude + co-latitude = 90 degrees.)" from: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html

To create the SPK kernel, you can get the inertial position vector by rotating the vector defined by the lat/long coordinates into an inertial frame, and then doing finite differences to calculate the velocity. This is how I did it in Python (using the SPICE python wrapper called spiceypy):

'''
Create SPK kernel from Lunar lat-lon coordinates
'''

from sys import path
path.append( '/home/alfonso/AWP/python_tools' )
import plotting_tools  as pt
import planetary_data  as pd
import spice_tools     as st
import numerical_tools as nt

import spiceypy as spice
import numpy as np

fn = '/mnt/c/Users/alfon/AWP/misc/'
r2d = np.pi / 180.0

if __name__ == '__main__':
spice.furnsh( '../../spice/lsk/naif0012.tls' )
spice.furnsh( '../../spice/pck/pck00010.tpc' )

lat     = 50.0 * r2d
lon     = 100.0 * r2d
r       = 1737.1
et0     = spice.str2et( '2021-05-17' )
etf     = spice.str2et( '2022-01-01' )
delta_t = 1000
ets     = np.arange( et0, etf, delta_t )
states  = np.zeros( ( len( ets ), 6 ) )
pos_iau_moon = spice.latrec( r, lon, lat )

for n in range( len( ets ) ):
states[ n, :3 ] = np.dot(
spice.pxform( 'IAU_MOON', 'J2000', ets[ n ] ),
pos_iau_moon )

for n in range( len( ets ) - 1 ):
states[ n, 3: ] = ( states[ n + 1, :3 ] - states[ n, :3 ] ) / delta_t

pos_np1 = np.dot(
spice.pxform( 'IAU_MOON', 'J2000', ets[ -1 ] + delta_t ),
pos_iau_moon )
states[ -1, 3: ] = ( pos_np1 - states[ -1, :3 ] ) / delta_t

handle   = spice.spkopn( 'lunar_point.bsp', 'SPK_file', 0 )
point_id = 301100
center   = 301
frame    = 'J2000'
degree   = 5
spice.spkw09( handle, point_id, center, frame,
ets[ 0 ], ets[ -1 ], '0', 5, len( ets ),
states.tolist(), ets.tolist() )
spice.spkcls( handle )

spice.furnsh( 'lunar_point.bsp' )
states_j2000    = st.calc_ephemeris( point_id, ets, 'J2000', 301 )
states_iau_moon = st.calc_ephemeris( point_id, ets, 'IAU_MOON', 301 )

latlons_j2000    = nt.inert2latlon( states_j2000[ :, :3 ], 'J2000', 'IAU_MOON', ets )
latlons_iau_moon = nt.bf2latlon( states_iau_moon[ :, :3 ] )

pt.plot_groundtracks(
[ latlons_j2000, latlons_iau_moon ],
{
'labels'         : [ 'J2000', 'IAU_MOON' ],
'plot_coastlines': False,
'filename'       : fn + 'lunar_point_groundtracks.png',
'dpi'            : 300
} )

pt.plot_orbits( [ states_j2000[ :, :3 ], states_iau_moon[ :, :3 ] ],
{
'labels'   : [ 'J2000', 'IAU_MOON' ],
'cb_cmap'  : 'gist_yarg',
'azimuth'  : 41.0,
'elevation': 30.0,
'axes_mag' : 1.0,
'filename' : fn + 'lunar_point_3d.png',
'dpi'      : 300
} )


In this case I used IAU_MOON but you could just replace that with MOON_ME and the ephemeris times that you are interested in.

The imports at the top (plotting_tools, planetary_data, spice_tools, numerical_tools) are all personal libraries, but you don't need them for the main part of the program (writing the BSP kernel), only after to get some visual confirmation that the BSP was written correctly. But for reference, here is the calc_ephemeris function in spice_tools (just a convenience wrapper to spkezr and spkgeo):

def calc_ephemeris( target, ets, frame, observer ):
'''
Wrapper for spkezr and spkgeo
'''

if type( target ) == str:
return np.array( spice.spkezr( target, ets, frame, 'NONE', observer )[ 0 ] )

else:
n_states = len( ets )
states   = np.zeros( ( n_states, 6 ) )
for n in range( n_states ):
states[ n ] = spice.spkgeo( target, ets[ n ], frame, observer )[ 0 ]
return states


I have part of the plotting_tools file posted here under src/python_tools: https://github.com/alfonsogonzalez/AWP

I show how to write a bunch of orbital mechanics software with Python in YouTube videos (if you just search "orbital mechanics with python" on YouTube you will find them, not sure if I'm allowed to post links).

And you're right, you could write this BSP kernel with the position w.r.t a body-fixed frame (IAU_MOON, MOON_ME), and I am just realizing that that would definitely be easier. At first, I did it in the inertial frame because calculating derivatives in the inertial frame are straight forward and you don't have to take into account the angular velocity of the frame since its inertial. But actually in this case the velocity of the point on the surface is just 0, so you're right I think you should do that instead. I just redid it that way and got what we would expect:

• Thank you very much. What are those python modules numerical_tools, spice_tools and planetary_data? One thing that really confuses me, eg. "MOON_ME" is a body fixed frame. I would think that it is enough to specify the frame centre with respect to this frame, not inertial frame (SPICE should take care of that), that is, without having to calculate velocities. Commented May 18, 2021 at 13:08
• I just updated the response to answer these questions. You're right about that you can write the SPK kernel w.r.t MOON_ME and its simpler to do so Commented May 18, 2021 at 14:07