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

      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'

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)?


1 Answer 1


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 ] )

        [ 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_radius': 1737.1,
        'cb_cmap'  : 'gist_yarg',
        'azimuth'  : 41.0,
        'elevation': 30.0,
        'axes_mag' : 1.0,
        'filename' : fn + 'lunar_point_3d.png',
        'dpi'      : 300
        } )

enter image description here

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 ] )

        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: enter image description here

enter image description here

  • 2
    $\begingroup$ 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. $\endgroup$
    – atapaka
    Commented May 18, 2021 at 13:08
  • 1
    $\begingroup$ 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 $\endgroup$ Commented May 18, 2021 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.