Basic Satellite 3-D East-North-Up Right-Handed Position Vector for Skyfield External Calculations

[jd7378]

First:

Skyfield is fantastic, thanks for the awesome product!

Second:

Wanted to share my learning in case it helps save time for anyone else who failed to pay attention to detail like I did. I take a TLE string, read into a satellite object, and then plot the altitude and azimuth, and compare that to a 3-D scatter plot in python. Originally I failed to use the topocentric reference frame, so I was messing things up with positions that referenced the wrong center. After setting that straight, I looked into the nuance of a North-East-Up left-handed coordinate system, as my vector calculations with reference frame x,y,z positions assumed they were taking placing in an East - North - Up right-handed frame. Lastly, don't forget radians for your azimuth, degrees for your elevation plot, or else you'll be in the wrong quadrant... Code below.

  from matplotlib import pyplot as plt
  
  %matplotlib tkagg 
  #pop out plots so they can be turned/manipulated for easier viewing
  
  #Establish reference timescales for platform positions based on TLEs
  from skyfield.api import load, wgs84
  ts = load.timescale() 
  ESTDelta = datetime.timedelta(hours=-5) #Eastern Standard Time (UTC-5)
  EST = datetime.timezone(offset=ESTDelta,name='EST') #Eastern Standard Time timezone
  now = datetime.datetime(year=2025,month=2,day=24,hour=10,minute=14,second=6,tzinfo=EST) #Datetime anchored to EST timezone, 24 Feb for repeatable charts
  datetimes = [now+datetime.timedelta(seconds=ii) for ii in np.arange(0,62.5,.5)] #one minute of data at 0.5 second intervals
  sky_times = ts.from_datetimes(datetimes)
  
  starlink = {'TLE_LINE0': '0 STARLINK-30504',
    'TLE_LINE1': '1 57980U 23151S   25054.83174474  .00002971  00000-0  12323-3 0  9998',
    'TLE_LINE2': '2 57980  42.9997 210.2000 0001938 266.6158  93.4472 15.25955393 79069'} #unused dict for easy reading
  
  #actual byte array from string for satellite load
  starlink_str = bytearray(str('1 57980U 23151S   25054.83174474  .00002971  00000-0  12323-3 0  9998\n2 57980  42.9997 210.2000 0001938 266.6158  93.4472 15.25955393 79069'),'utf-8') #TLE for Starlink 30504, 24 Feb 2025
  
  from io import BytesIO
  f = BytesIO(starlink_str)
  
  from skyfield.iokit import parse_tle_file
  satellite = list(parse_tle_file(f, ts))
  
  #Define Earth Observation Point (CIS)
  cis = wgs84.latlon(43.086091468700196,-77.67769150880549,0.0) #CIS Building, with receiver located at wgs84 sea level (not accurate for real collect...)
  
  def ENU_NEU_Swap(vector):
      """
      Swap 1st and 2nd Component of a vector to modify coordinate system from ENU to NEU, or NEU to ENU.
      
      Parameters:
      Vector (array-like): Input vector in ENU coordinate system
      
      Returns:
      Vec[float]: Modified vector in NEU coordinate system
      """
      vector[0],vector[1] = vector[1],vector[0]  # Swap East and North components
      
      return vector
  
  print((satellite[0]-cis).at(sky_times[15]).frame_xyz(cis).km)
  
  temp_ax = plt.figure().add_subplot(projection='3d')
  temp_ax.scatter(*[0,0,0], label='CIS')
  sat_position = (satellite[0]-cis).at(sky_times[15]).frame_xyz(cis).km #get satellite position in km
  sat_position_ENU = ENU_NEU_Swap(sat_position) #convert to ENU coordinates
  temp_ax.scatter(*sat_position_ENU,label=satellite[0].name)
  temp_ax.scatter(*[100,100,0]) #Placeholder points to make visual cleaner
  temp_ax.scatter(*[-100,100,0]) #Placeholder points to make visual cleaner
  temp_ax.scatter(*[100,-100,0]) #Placeholder points to make visual cleaner
  temp_ax.scatter(*[-100,-100,0]) #Placeholder points to make visual cleaner
  temp_ax.quiver(*np.zeros(3,), *[0,100,0], color='blue', label='North Vector')
  temp_ax.quiver(*np.zeros(3,), *[100,0,0], color='red', label='East Vector')
  temp_ax.set_aspect('equal')
  temp_ax.legend()
  temp_ax.set(xlabel='X (km)', ylabel='Y (km)', zlabel='Z (km)')
  temp_ax.view_init(elev=90, azim=-90, roll=0)
  temp_ax.set(title=f'Satellite 3-D Position with ENU Swap')
  
  temp_ax = plt.figure().add_subplot(projection='3d')
  temp_ax.scatter(*[0,0,0], label='CIS')
  sat_position = (satellite[0]-cis).at(sky_times[15]).frame_xyz(cis).km #get satellite position in km
  temp_ax.scatter(*sat_position,label=satellite[0].name)
  temp_ax.scatter(*[100,100,0]) #Placeholder points to make visual cleaner
  temp_ax.scatter(*[-100,100,0]) #Placeholder points to make visual cleaner
  temp_ax.scatter(*[100,-100,0]) #Placeholder points to make visual cleaner
  temp_ax.scatter(*[-100,-100,0]) #Placeholder points to make visual cleaner
  temp_ax.quiver(*np.zeros(3,), *[0,100,0], color='blue', label='North Vector')
  temp_ax.quiver(*np.zeros(3,), *[100,0,0], color='red', label='East Vector')
  temp_ax.set_aspect('equal')
  temp_ax.legend()
  temp_ax.set(xlabel='X (km)', ylabel='Y (km)', zlabel='Z (km)')
  temp_ax.view_init(elev=90, azim=-90, roll=0)
  temp_ax.set(title=f'Satellite 3-D Position without ENU Swap')
  
  
  fig, temp_ang = plt.subplots(subplot_kw={'projection': 'polar'},layout='constrained')
  temp_ang.plot((satellite[0]-cis).at(sky_times[15]).altaz()[1].radians, (satellite[0]-cis).at(sky_times[15]).altaz()[0].degrees, '*') 
  temp_ang.set_theta_zero_location('N')
  temp_ang.set_theta_direction(-1)
  temp_ang.set_rlim([90,40])
  plt.rcParams['axes.titley'] = 1.0    # y is in axes-relative coordinates.
  plt.rcParams['axes.titlepad'] = 20
  temp_ang.set(title=f'Starlink Satellite Horizonal Reference')