Rotation_at(): Produce rotation matrix from ICRS to topocentric frame at Geographic Position

[concord647]

Hello again,

I need a way to generate the rotation matrix from ICRS to a (topocentric) coordinate system (that may or may not be defined as shown in the drawing) defined at a greographic location (i.e. a Ground Station) for specific times.

I am looking at Skyfield's rotation_at()-method, but I don't fully understand it. Say a produce a geographic position like this:

GS_Lat = 34.22475000  # North
GS_Lon = 118.05591667  # West
GS_Alt = 1728  # meter

GS = wgs84.latlon(GS_Lat * N, GS_Lon * W, elevation_m=GS_Alt)

I can then produce a rotation via rotation_at for a specific time:

ts = load.timescale()
timenow = ts.now()

# Geocentric position in ICRS at timenow
position_now = GS.at(timenow)

# Rotation ICRS --> what frame??
rotation_now = GS.rotation_at(timenow)

This produces some rotation matrix from ICRS --> GS. What I am slightly confused about is the geographic coordinate system. How is this defined, i.e. what exactly is this rotation matrix rotating the ICRS frame to?

Thanks a lot!

[brandon-rhodes]

Here's how the documentation defines that reference frame:

https://rhodesmill.org/skyfield/coordinates.html#horizonal-coordinates

Is there a way the title of that section in the docs could be tweaked or expanded to have made it easier for you to find? Should "altitude" or "azimuth" be added to the text of the title, for example?

[brandon-rhodes]

When I made my diagrams, that was a symptom of my applying the matrix backwards. Once I used the transpose, they all lined up north-east-up. Could you check whether you are using the .T of the matrix?

[concord647]

Yes, sorry I should have stated this: the new plots above (with the colorful earth surface) are with the transposed rotation matrix of the rotation_at()-output, using your T-function:

r1 = R.from_matrix(T(GS.rotation_at(time_vector[int(i)])))  # Transposed rotation matrice --> GS --> ICRS

[brandon-rhodes]

If I had more time at the moment I would clean this up a bit, but here's how I did the rotation in my edited version of your script. It looks like I used my mxv() matrix-times-vector call, and used plain-old .T because (whoops) I forgot about my T() function. See if the same maneuver works for you — this script gives me what look like upright axes; does it work for you too on your system?

from skyfield.api import load, wgs84, N, S, E, W
from skyfield import api, units, framelib
from skyfield.functions import T, mxv

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial.transform import Rotation as R

def plot_rotated_axes(ax, r, rot, name=None, offset=(0, 0, 0), scale=2e6):
    print('k')
    colors = ("#FF6666", "#005533", "#1199EE")  # Colorblind-safe RGB
    labels = ('x', 'y', 'z')
    loc = np.array([offset, offset])

    for i, (axis, c, ax_label) in enumerate(zip((ax.xaxis, ax.yaxis, ax.zaxis), colors, labels)):
        axlabel = ax_label
        axis.set_label_text(axlabel)
        axis.label.set_color(c)
        axis.line.set_color(c)
        axis.set_tick_params(colors=c)
        if rot is None:
            line = np.zeros((2, 3))
            line[1, i] = scale
            line_rot = r.apply(line)
        else:
            v = np.array((0.0, 0.0, 0.0))
            v[i] = scale
            v = mxv(rot.T, v)
            line_rot = np.zeros((2, 3))
            line_rot[1] = v
        line_plot = line_rot + loc
        ax.plot(line_plot[:, 0], line_plot[:, 1], line_plot[:, 2], c)
        text_loc = line_rot[1]*1.2
        text_loc_rot = r.apply(text_loc)
        text_plot = text_loc_rot + loc[0]
        ax.text(*text_plot, axlabel.upper(), color=c,
                va="center", ha="center")
    # ax.text(*offset, name, color="k", va="center", ha="center",
            # bbox={"fc": "w", "alpha": 0.8, "boxstyle": "circle"})

## DEFINE GEODETIC OGS COORDINATE SYSTEM
GS_Lat = 0
GS_Lon = 0
GS_Alt = 0

GS = wgs84.latlon(GS_Lat * N, GS_Lon * W, elevation_m=GS_Alt)

ts = load.timescale()
timenow = ts.now()

# Geocentric position in ICRS at timenow
position_now = GS.at(timenow)

# Rotation ICRS --> what frame??
rotation_now = GS.rotation_at(timenow)

print(position_now.position.m)

## Test: Get ground station position in ICRS over full day and plot in 3D
# define start and end time in julian datum
time_start = ts.tt(year=2023, month=10, day=24)
time_end = ts.tt(year=2023, month=10, day=25)
steps = 24

time_vector = ts.linspace(time_start, time_end, num=steps)

# get position of ground station at time_vector
position_vector_fullday = GS.at(time_vector)

# plot
fig = plt.figure()
ax = fig.add_subplot(projection='3d', proj_type="ortho")

# plot ground station positions
ax.scatter(position_vector_fullday.position.m[0], position_vector_fullday.position.m[1], position_vector_fullday.position.m[2], color="k", s=20)

# add sphere to plot
earth_radius = 6371*1000+GS_Alt

# Make data for sphere
u, v = np.mgrid[0:2*np.pi:50j, 0:np.pi:50j]
x = earth_radius*np.cos(u)*np.sin(v)
y = earth_radius*np.sin(u)*np.sin(v)
z = earth_radius*np.cos(v)

# Plot earth surface
ax.plot_surface(x, y, z, rstride=1, cstride=1, color='c', alpha=0.3, linewidth=0)

# plot reference axis
r_none = R.from_matrix(np.identity(3))
plot_rotated_axes(ax, r_none, None, name="r1", offset=[0,0,0])
#plt.show()

# plot horizonal axis at surface points
plot_fraction = 6
plot_indices = np.linspace(0, steps, int(steps/plot_fraction)+1)
print(plot_indices)

for i in np.nditer(plot_indices[0:len(plot_indices)-1]):    
    # r1 = R.from_matrix(np.identity(3))  # translation of CS without rotation
    # r1 = R.from_euler("ZYX", [90, 0, 0], degrees=True)  # single rotation about Z (90°)
    
    # r1 = R.from_matrix(T(GS.rotation_at(time_vector[int(i)])))  # Transposed rotation matrice --> GS --> ICRS ??
    
    rot = GS.rotation_at(time_vector[int(i)])
    r1 = R.from_matrix(rot)  # rotation: ICRS --> GS-Frame (Horizonal)

    # euler angles
    euer_angles = r1.as_euler('ZYX', degrees=True)
    print(f'Euler-Angles (ZYX): {euer_angles}')

    # check orthogonality
    print(f'Det of Rotation: {np.linalg.det(GS.rotation_at(time_vector[int(i)]))}')
    
    plot_rotated_axes(ax, r1, rot, name="r1", offset=position_vector_fullday[int(i)].position.m)
    
plt.show()

[concord647]

Thanks Brandon, this does indeed solve the issue:

It took me a while to understand why this produces a different outcome than using Scipy.spatial.transform.Rotation, but I think I got it now: I hadn't realised that the ouput of GS.rotation_at() has a negative determinant, incorporating the "flip" from right-handed to left-handed. R.from_matrix() according to the scipy-docs checks for det=1 and "If the input is not proper orthogonal, an approximation is created using the method described in [2]." This approximation changes the transformation in some way.

If I "unflip" the output of GS.rotation_at() with [[1,0,0],[0,-1,0],[0,0,1]], then use the inverse of that with R.apply() for actual rotating , I get a right-handed topocentric coordinate system ("North-West-Up" I guess):

I completely overlooked that your Horizonal definition is left-handed! Is there a particular reason for that?

Thanks a lot for your support!

[brandon-rhodes]

I completely overlooked that your Horizonal definition is left-handed! Is there a particular reason for that?

Yes, it's because the associated angles — azimuth and altitude — measure "longitude" clockwise around the unit circle instead of counter-clockwise. I have seen systems that instead of returning an altitude angle (also called elevation) instead return a "nadir angle" or somesuch, that's negative for objects above the horizon. But as so many Skyfield users are amateur astronomers, I wanted to go with the standard convention of positive "latitude" being up towards the zenith.

I'm glad we now have your axes correctly oriented! I had never heard of Scipy.spatial.transform.Rotation before, but now I'll know about a possible failure condition if I ever see someone using it again on a problem like this.