1
$\begingroup$

Original question below, but my question is now about how can I supply lat/lon coordinates for other locations other than the one in the example below?

It is possible to use NASA's SPICE toolkit to produce Alt/Az coordinates? I haven't been able to find any examples.

Update:2 The code below doesn't appear to be along the right lines. I've posted an answer with the "most" correct code I could come up with.

Update:

I think I partly figured it out. Apparently most people call it Elevation and Azimuth, not AltAz, which is why I couldn't find anything. I found a SpiceyPy example which generates the Az, El for DSS-13. But I'm at a loss as to how to specify the Lat,Lon for general locations. Looking through the kernel file, the closest thing I see that could be a lat/lon for DSS-13 are the "TOPO_ANGLES" which are -243.2055404763616,-54.7528357325366, which is a pretty odd format if its supposed to be lat,lon. I was unable to find the real lat/lon for DSS-13, but I was able to find it for DSS-14 which appears in the same file (FYI those coords are 116.89W, 35.426666666666N). I verified the results with Horizons and they match exactly.

Below I provide a stripped down version of the example to only print out the el/az coordinates. It requires a file "erotat.tm" which I also provide below. I was able to filed all of the required kernel files by searching for the name (apparently some of them are quite old). I was not able to find de414_2000_2020.bsp, but I was able to use de414.bsp from here. You can install SpiceyPy (assuming you already have Python and PIP) by typing "pip install spiceypy".

import spiceypy

METAKR = 'erotat.tm'
spiceypy.furnsh( METAKR )
et = spiceypy.str2et('2007 JAN 1 00:00:00')
[topov, ltime] = spiceypy.spkpos( 'moon', et, 'DSS-14_TOPO','lt+s', 'DSS-14')
[r, lon, lat] = spiceypy.reclat( topov )

az = -lon
if  az < 0.0:
    az += spiceypy.twopi()
el = lat

print( 'Moon Az (deg):        {0:15.6f}\n'
       'Moon El (deg):        {1:15.6f}\n'.format(
           az * spiceypy.dpr(),
           el * spiceypy.dpr() )  )

spiceypy.unload( METAKR )

Contents of emrotat.tm

KPL/MK

Meta-kernel for the "Earth Rotation" task
in the Binary PCK Hands On Lesson.

The names and contents of the kernels referenced by this
meta-kernel are as follows:

File name                       Contents
------------------------------  ---------------------------------
naif0008.tls                    Generic LSK
de414_2000_2020.bsp             Solar System Ephemeris
earthstns_itrf93_050714.bsp     DSN station Ephemeris
earth_topo_050714.tf            Earth topocentric FK
pck00008.tpc                    NAIF text PCK
earth_000101_070725_070503.bpc  Earth binary PCK


\begindata

KERNELS_TO_LOAD = ( 'naif0008.tls'
                    'de414.bsp'
                    'earthstns_itrf93_050714.bsp'
                    'earth_topo_050714.tf'
                    'pck00008.tpc'
                    'earth_000101_070725_070503.bpc' )

\begintext
$\endgroup$
  • 1
    $\begingroup$ I'm not a spice user so I can't answer, but I'm pretty sure the answer will turn out to be "yes". I use the Horizons website to calculate altitude and azimuth of objects from the surface of other solar system objects, here's one example so if Horizons can do it, then Spice can. $\endgroup$ – uhoh Aug 20 at 0:24
  • 1
    $\begingroup$ My thoughts exactly, I'd think it'd be the #1 function done with it, but I've searched a lot and can't find any examples or even discussion about it. $\endgroup$ – Greg Miller Aug 20 at 1:59
  • $\begingroup$ If you think you have even a partial answer (it looks pretty complete to me) that would be potentially beneficial to future readers, you should post it as an answer post, rather than an update to the question. That way people will see that an answer exists and come and read it. As a side benefit you pick up a few reputation points as well. It is always okay to answer your own question, and if nothing happens in a few days, accept it as well. $\endgroup$ – uhoh Aug 20 at 4:02
  • 1
    $\begingroup$ I almost did, but it's only been about a day, hopefully someone will chime in with the whole thing. Otherwise I'll post it once I have the whole thing figured out. $\endgroup$ – Greg Miller Aug 20 at 4:10
2
$\begingroup$

I kinda found a solution, unfortunately it's a few hundredths of a degree off from those specified by Horizons. And it requires you to define a custom kernel for the observer's reference frame, which isn't bad as long as your doing a lot of observations from one location. The code below is in Python and requires spiceypy ("pip install spiceypy").

import spiceypy
import math
#By Greg Miller (gmiller@gregmiller.net) 2019
#Released as public domain

spiceypy.furnsh('erotat.tm')

#Lat,Lon position of observer e.g. Louisville, KY
#Needs to match definition of "MYTOPO" in custom kernel mytopo.tf
lat=math.radians(38.2464000)
lon=math.radians(274.236400)

#convert lat,lon to xyz
obspos = spiceypy.georec(lon, lat, 0.0, 6378.1366, 1.0/298.25642)

#time of observation in UT
et = spiceypy.str2et('2007 JAN 1 00:00:00')

#Get the xyz coordinates of the moon relative to observer in MYTOPO frame
moon, lt = spiceypy.spkcpo('MOON', et, 'MYTOPO', 'OBSERVER', 'LT+S', obspos, 'EARTH', 'ITRF93')

#reclat function only wants a length 3 array
temp=moon[0:3]

#convert to polar coordinates
[r, lon, lat] = spiceypy.reclat( temp )

#Convert polar coords to Az, El
az = -lon

if  az < 0.0:
    az += spiceypy.twopi()

el = lat

print( 'Moon Az (deg):        {0:15.6f}\n'
       'Moon El (deg):        {1:15.6f}\n'.format(
           az * spiceypy.dpr(),
           el * spiceypy.dpr() )  )


spiceypy.unload('erotat.tm')

This is the contents of the erotate.tm kernel spec file. You must also download all of the kernel files individually, you can generally find them by searching for their names. Note that the last one "mytopo.tf" is one you must supply, an example is below.

\begindata

KERNELS_TO_LOAD = ( 'naif0008.tls'
                    'de414.bsp'
                    'earthstns_itrf93_050714.bsp'
                    'earth_topo_050714.tf'
                    'pck00008.tpc'
                    'earth_000101_070725_070503.bpc',
                    'mytopo.tf' )

\begintext

Below is an example mytopo.tm file. The lon,lat go in the TKFRAME_1234567_ANGLES and is the only field you'll need to edit. Note the odd format for lat/lon compared to how it's specified in the Python code above. The lon is just the negative, the lat is the negative of 90-lat.

\begindata

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

  TKFRAME_1234567_SPEC     = 'ANGLES'
  TKFRAME_1234567_RELATIVE = 'IAU_EARTH'
  TKFRAME_1234567_ANGLES   = ( -274.236400, -51.7536, 180 )
  TKFRAME_1234567_AXES     = (       3,        2,   3 )
  TKFRAME_1234567_UNITS    = 'DEGREES'

\begintext
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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