Is it posible to convert JPL Horizons Vectors to ECEF?

Question

I'm needing to accurately find the sun and moon position in ECEF and I'm wanting to make use of JPL's Horisons VECTORS to do this.

What data should I export from the web interface and how should I convert it?

---

Background

I need the accurate positions to find 3d angles of incoming light onto slopped surfaces and also if the view of a point on that surface of the sun or moon is obstructed by another location on that surface that i have mapped in ECEF.

There are around 500 surfaces about 20sqkm each derived from LiDAR sources located at various parts of the earth.

---

What worked for me

I ended up following @Alfonso Gonzalez's advice and getting the position directly from SPICE, and I'm able to get the same results as the Horizons OBSERVER data but doing the math in the ITRF93/ECEF reference frame.

This is the python code I'm using:

import spiceypy as spice

spice.furnsh('earth_000101_210617_210326.bpc')
spice.furnsh('naif0012.tls')
spice.furnsh('de432s.bsp')

et = spice.str2et('2021-03-23T00:00:00')

sun_naif_id = spice.bods2c('SUN')
earth_naif_id = spice.bods2c('EARTH')

sun, lt = spice.spkezp(sun_naif_id, et, 'ITRF93', 'LT+S', earth_naif_id)
sun *= 1000

print("Sun Position at 2021-03-23 00:00 UTC in ECEF\n  ", sun)

Outputs:

Sun Position at 2021-03-23 00:00 UTC in ECEF
   [-1.49005330e+11 -4.33523011e+09  2.67037253e+09]

As test I tried calculating the sun elevation angle for Adelaide, Australia at 2020-03-23 00:00:00 UTC, and got 36.060167 and Horizons OBSERVER gives 36.0602 for the same time and location.

import numpy as np

adelaide = np.array([-3927127.15169169, 3462288.85408112, -3630705.11577511])
adelaide_z = np.array([-4542168.38388766, 4004530.12124128, -4203154.84360054])

ba = adelaide_z-adelaide
bc = sun-adelaide
cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
z_angle = np.arccos(cosine_angle)

print("Sun Elevation for Adelaide 2021-03-23 00:00 UTC\n  ", 90 - np.degrees(z_angle))

Outputs:

Sun Elevation for Adelaide 2021-03-23 00:00 UTC
   36.06016701724842

Answer 1

You could do all this analysis using SPICE, including getting ephemeris data from SPICE kernels instead of Horizons. I'm not sure which programming language you're using, but you can get the toolkit directly in C, Fortran, MATLAB, and Java (https://naif.jpl.nasa.gov/naif/toolkit.html) or there's a 3rd party Python wrapper that I use all the time called SpiceyPy (https://spiceypy.readthedocs.io/en/master/documentation.html)

First, if you want to get the ephemeris from SPICE, you'll need an ephemeris kernel, which can be found here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/ de432s.bsp will work for what you're doing. And then you'll need a leapseconds kernel, which can be found here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/ naif0012.tls will work.

In order to convert from inertial to ECEF vectors, there are several options for Earth orientation kernels, but I like to use earth_200101_990628_predict.bpc most of the time, which can be found here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/

Once you have all these, you can use the SPICE toolkit to get the ECEF vectors. Here is how I did it in Python for the groundtrack plot:

If you want to use the data you got from Horizons and convert it to ECEF, you can use SPICE's pxform function to obtain the rotation matrices at each time step. That would look something like this:

If you want to do this in python, you can install spiceypy using pip:

$ python3 -m pip install spiceypy

The conversion to latitude / longitude coordinates is out of the scope of the question, but if you're interested, that bf2latlong function is using SPICE's reclat function, which converts from rectangular to latitudinal coordinates. https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/reclat_c.html

Also, SPICE has functions for calculating illumination angles and local solar time, which you could also use for your analysis: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ilumin_c.html

Answer 2

This is a simple answer that is not super-accurate but it will get you started.

Hopefully someone will chime in with corrections/refinements!

---

If we work in ecliptic coordinates then for a spherical Earth of radius $r_0$ the position of a location is just

\begin{align} x & = r_0 \cos\phi \cos(\lambda + \theta)\\ y & = r_0 \cos\phi \sin(\lambda + \theta)\\ z & = r_0 \sin\phi \\ \end{align}

where $\phi$ and $\lambda$ are latitude and longitude, $\theta$ is the rotational angle of the Earth.

Presumably $x, y, z$ will be arrays in your program, of the same size and dimensionality as $\phi$ and $\lambda$.

If you want to use an oblate Earth and include height above the WGS84 geoid (i.e. altitude) the see this answer.

One approximation here is that the precession of Earth's axis is not taken into account. Below you'll see that the epoch for the Earth's equator is J2000 and it's moved a tiny bit since then.

From Earth rotation angle:

$$\theta(t_U) = 2 \pi (0.779\,057\,273\,2640+1.002\,737\,811\,911\,354\,48t_{U})$$

where $t_U$ is the Julian UT1 date − 2451545.0.

For this approximation you can use any Julian date converter your programming environment provides. For example the Julian date of 27-Mar-2021 12:00 UTC is about 2459301.00. Note that JD crosses zero integer value at noon, not midnight, so as not to confuse poor astronomers in the middle of the night.

Please be sure to ask Horizons to give you the positions of the Sun, Moon, and Earth in

reference plane: Earth mean equator and equinox of reference epoch
reference system: ICRF/J2000.0

per this answer. That looks like this:

Then subtract the Sun's position from your $x, y, z$ array of points on the Earth to get vectors at each point that point towards the Sun, and similar for the Moon. You should see as much as a +/- 1° shift in those vectors for the Moon at opposite sides of the Earth because the Moon is so close.

Then once you have the vectors from the faces to the Sun or to the Moon, you can normalize them.