An implementation question about altaz

[karatemir]

Hi @brandon-rhodes

I was testing out the altaz functionality and have some trouble with the implementation. I hope you can clarify some of my confusion. As a test scenario let us consider the following:

import numpy as np
from skyfield.api import EarthSatellite, load, wgs84
from datetime import datetime, timedelta

tle = """\
1 39444U 13066AE  22267.86615744  .00003125  00000+0  37269-3 0  9990
2 39444  97.6351 237.3665 0057839  66.8984 293.8312 14.83771785476266"""

sat = EarthSatellite(*tle.splitlines())
t0 = sat.epoch + timedelta(minutes=15)

print(wgs84.subpoint_of(sat.at(t0)))
# prints WGS84 latitude +55.5663 N longitude -93.0798 E elevation 0.0 m

loc = wgs84.latlon(55, -90)
alt, az, dist = (sat - loc).at(t0).altaz()

print(f'elevation: {alt.degrees} , azimuth: {az.degrees}, slant range: {dist.km}')
# prints elevation: 68.95630788820152 , azimuth: 289.1205055413545, slant range: 624.9926563478068

I then tried reproducing the same numbers using basic vector geometry. The idea is we have two vectors: Geocentric vector for the Earth station loc and another Geocentric vector for sat. Taking the difference, then the norm should ideally give us the slant range:

v = loc.at(t0).position.km
w = (sat - loc).at(t0).position.km
print(f'slant range: {np.linalg.norm(w)}')
# slant range: 624.9926563478072

So altaz and vector geometry calculations are basically equal up to machine precision, so far so good. To calculate the elevation angle we can calculate the dot product between w (difference vector) and v (location vector).

el = np.rad2deg (np.arccos ( np.dot(w,v) / np.linalg.norm(w) / np.linalg.norm(v) ) )
el = 90 - el # elevation angle is defined from the horizon
print(f'elevation angle: {el}')
# prints elevation angle: 68.89634982362499

So we have ~0.06 degree discrepancy in elevation angle. To compute the azimuth angle I follow the discussion here.

d = v / np.linalg.norm(v)
z = np.array([0,0,1])
e = np.cross(d,z)
e = e / np.linalg.norm(e) # normalize east vector
n = np.cross(d,e) # north pointing vector
p = np.cross(w,d) # this vector will be at 90 degrees offset to the projection vector to the tangential plane
p = p / np.linalg.norm(p)
minor = np.linalg.det( np.stack( (n[-2:], p[-2:]) ) ) # for -pi to pi coverage
sign = 1 if minor == 0 else np.sign(minor)
dot_p = np.clip(np.dot(n, p), -1.0, 1.0)
az = np.rad2deg(sign * np.arccos (dot_p ))
print(f'azimuth: {(az + 90) %360}')
# azimuth: 289.7260847243732

So in azimuth calculation we have ~0.6 degrees discrepancy. As a sanity check I also used pyproj to get the geodetic azimuth angle.

from pyproj import Geod
g = Geod(ellps='WGS84')
s = wgs84.subpoint_of(sat.at(t0))
az, _, _ = g.inv(loc.longitude.degrees, loc.latitude.degrees, s.longitude.degrees, s.latitude.degrees)
print(f'pyproj azimuth: {az % 360}')
# prints pyproj azimuth: 289.1174382871312

This value is within 0.003 degrees of SkyField azimuth value so much better than my computation.

Digging through the code, I see that altaz is computed via here and here. Basically after applying a rotation, the spherical coordinates are reported as alt, az, dist values.

Can you please explain what is incorrect with my approach? Am I missing something?

[brandon-rhodes]

Can you please explain what is incorrect with my approach? Am I missing something?

The vector from the geocenter to a point on the Earth's surface is not coincident with the vector to the zenith at that point, because the Earth is not a sphere but an oblate spheroid. (Unless the point is at the equator or the poles.)

https://en.wikipedia.org/wiki/Geodetic_coordinates#Geodetic_vs._geocentric_coordinates

I think that might explain why your angles are coming out a bit off?

[karatemir]

@brandon-rhodes thanks for prompt response. Is there a built-in way to get the geodetic normal vector (zenith pointing) ?

[brandon-rhodes]

Is there a built-in way to get the geodetic normal vector (zenith pointing) ?

Try asking the position for its .frame_xyz(loc) and see if a vector comes back out.

[karatemir]

Just wanted to update:

Following your advice I was able to reproduce altaz values through geometric vector operations.

To get the normal I used

loc2 = wgs84.latlon(loc.latitude.degrees, loc.longitude.degrees, elevation_m=1000)
v = (loc2 - loc).at(t0).position.km

I tried loc.at(t0).frame_xyz(loc) as a more principled alternative but I couldn't make sense out of it.

I was also using [0,0,1] for unit vector pointing towards the spin axis of Earth. Later on I figured out that GCRF z is not coincident with spin axis. To get spin axis I use: z = wgs84.latlon(90,0).at(t0).position.km. With these changes in place things worked great. Thanks once again for your advice.

As a final step I also tried out from_altaz to go back to the satellite position from Earth Station. At first I wasn't able to understand what was going on but this stackoverflow question that you answered was very helpful. Maybe it would be useful to have some of that information in the documentation.

[brandon-rhodes]

Following your advice I was able to reproduce altaz values through geometric vector operations.

Good, I'm glad you were able to produce them! Have you checked out the rotation matrix itself, by the way? My intuition is that the three unit vectors are there inside the rotation matrix. Try either:

R = loc.rotation_at(t0)

—or else:

R = loc.rotation_at(t0).T

—and then examine the vectors:

print(R[0])
print(R[1])
print(R[2])

Do any of those vectors look like the unit vectors you're looking for?

[karatemir]

v = (loc2 - loc).at(t0).position.km
>>> v
array([0.49580111, 0.07782549, 0.86494188])

>>> loc.rotation_at(t0)
array([[-8.54382680e-01, -1.34739684e-01,  5.01871950e-01],
       [-1.55600425e-01,  9.87820030e-01,  3.11335412e-04],
       [ 4.95801114e-01,  7.78254890e-02,  8.64941876e-01]])

looks like R[2] does the job! could be exposed to user as .normal maybe

[brandon-rhodes]

could be exposed to user as .normal maybe

Maybe—but for other reference frames the name would need to be different, like .north_pole or .ecliptic_pole. Maybe instead the documentation could show examples of unpacking the rotation matrix into its three vectors? Just making things up, not yet checking the order they would really arrive in, it might look like:

north, east, up = loc.rotation_at(t0)
north_pole, zero_longitude, ninety_longitude = itrs.rotation_at(t0)

And so forth.

[karatemir]

that example is very useful! I confirmed that my north and east calculations match rotation_at unpacking (NEU as you suggested). Would be great if documentation can include a table corresponding to what each row corresponds to for different frames.