import csv
import matplotlib.pyplot as mplot
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import copy
import urllib.request as urls
import ssl
from astropy.time import Time as ast_time
[docs]class ER2_Flight_Data(object):
def __init__(self, source: str, verbose: bool=False, bore_offset_horz: float=0, bore_offset_vert: float=-3.23,
bore_offset_top: float=2, vert_fov: float=4, num_los: int=10, time_range=[]) -> None:
"""
Constructs the ER2_Flight_Data object. This object retrieves the flight data from
https://asp-archive.arc.nasa.gov and produces the attitude solution of the aircraft
and lines of sight of the instrument.
This is the only core function that accepts direct input form the user. Once instantiated all properties will
be populated automatically.
:param source: (string) The text specifying the source of the flight data. This is one of the following inputs:
- The url link to a specific flight, example:
"https://asp-archive.arc.nasa.gov/N809NA/FY2017/2017-04-20/IWG1.21Apr2017-0459"
- The flight simulation to use instead of a real one, example:
"SimOne" will load Flight_Simulation_One.csv,
"SimTwo" will load Flight_Simulation_Two.csv, etc
- A one line string which has flight information for only a single point in time.
This is used for real time input, example:
"IWG1,2017-07-19T18:02:15.000,34.1414141414,-118.070707071,22000.0, ......
:param verbose: (bool, optional, default=False) If true will print out statements letting the user know what
is being done.
:param bore_offset_horz: (float) The angle in degrees between the bore sight of the instrument and the forward
vector of the aircraft (the nose) in the lateral plane of the aircraft. Defined as
positive pointing to the right when facing forward.
:param bore_offset_vert: (float) The angle in degrees between the bore sight of the instrument and the forward
vector of the aircraft (the nose) in the normal plane of the aircraft. Defined as
positive pointing upward.
:param bore_offset_top: (float) The angle in degrees between the bore sight of the instrument and the upper most
field of view of the instrument (ideally this is half of the total fov).
Defined as always positive.
:param vert_fov: (float) The angle in degrees of the total field of view.
:param num_los: (float) The number of lines of sight to produce (ideally equal to the number of rows on
the detector).
:param time_range: (list, optional, default=[]) A list which defines the time, or range or times, to
calculate the geocentric attitude, position, and lines of sight for. This is done to
allow the user to save computation time. If the list is empty (default) data for all
times will be calculated. If the list is None no calculations will be done, only the raw
ER2 data will be read in. If the list is of length one only the closest time to that
given time will be calculated. If the list is of length two then these two times will
define the range of times to look at. The time into the list should be in hours in 24
hours UTC.
Author: Daniel Letros
Date Revised: 2017-05-08
"""
self.__verbose__ = verbose
# ============= Instantiating Properties =============
#
# ----------Raw flight data from url----------
# All of these properties are numpy arrays of equal length.
self.TimeStamp = []
self.Time_Hrs = []
self.mjd = []
self.Latitude = [] # [+- 90 deg]
self.Longitude = [] # [+- 180 deg]
self.GPS_Altitude_MSL = [] # [m]
self.GPS_Altitude = [] # [m]
self.Pressure_Altitude = []
self.RADAR_Altitude = [] # [ft]
self.Ground_Speed = [] # [m/s]
self.True_Air_Speed = [] # [m/s]
self.Indicated_Air_Speed = [] # [m/s]
self.Mach_Number = []
self.Vertical_Speed = [] # [m/s]
self.True_Heading = [] # [+- 180 deg] Positive CW from north, how nose is pointed from true north
self.Track_Angle = [] # [+- 180 deg] Positive CW from north
self.Drift_Angle = [] # [+- 180 deg] Positive CW from north
self.Pitch_Angle = [] # [+- 180 deg] Positive Up, nose up/down from level
self.Roll_Angle = [] # [+- 180 deg] Positive Right wing down, rotation left/right from level
self.Slip_Angle = [] # [deg]
self.Attack_Angle = [] # [deg]
self.Static_Air_Temp = [] # [C]
self.Dew_Point = [] # [C]
self.Total_Air_Temp = [] # [C]
self.Static_Pressure = [] # [mbar]
self.Dynamic_Pressure = [] # [mbar]
self.Cabin_Pressure = [] # [mbar]
self.Wind_Speed = [] # [0-132 m/s]
self.Wind_Direction = []
self.Vert_Wind_Speed = [] # [m/s]
self.Solar_Zenith_Angle = [] # [deg]
self.Aircraft_Sun_Elevation = [] # [deg]
self.Sun_Azimuth = [] # [deg]
self.Aircraft_Sun_Azimuth = [] # [deg]
# ----------Calculated Data----------
# All of these properties are numpy arrays
# index: [time_range, [x, y, z]]
#
# Note that internally the code populates these as list of lists.
# index: [time_range][x, y, z]
# These are turned into numpy arrays only as a final step in the function
# self.__convert_lists_to_ndarrays__()
self.aircraft_position_geocentric_xyz = [] # [m]
# Vectors which are local to aircraft position and with respect to Earth
self.aircraft_local_up_geocentric_xyz = [] # normalized vector
self.aircraft_local_down_geocentric_xyz = [] # normalized vector
self.aircraft_local_north_geocentric_xyz = [] # normalized vector
self.aircraft_local_south_geocentric_xyz = [] # normalized vector
self.aircraft_local_east_geocentric_xyz = [] # normalized vector
self.aircraft_local_west_geocentric_xyz = [] # normalized vector
# Vectors which are local to the aircraft position and with respect to its attitude
self.aircraft_attitude_up_geocentric_xyz = [] # normalized vector
self.aircraft_attitude_down_geocentric_xyz = [] # normalized vector
self.aircraft_attitude_forward_geocentric_xyz = [] # normalized vector
self.aircraft_attitude_backward_geocentric_xyz = [] # normalized vector
self.aircraft_attitude_left_geocentric_xyz = [] # normalized vector
self.aircraft_attitude_right_geocentric_xyz = [] # normalized vector
# ----------Instrument LOS Info----------
self.bore_offset_top = bore_offset_top # [deg] Always positive
self.bore_offset_horz = bore_offset_horz # [+- 180 deg] Positive right from nose of plane
self.bore_offset_vert = bore_offset_vert # [+- 180 deg] Positive up from nose of plane
self.vert_fov = vert_fov # [deg]
self.num_los = num_los
# self.instrument_los is a ndarray [time_range, bottom:top fov, [x,y,z]]
#
# Note that internally the code populates this as a list [time_range][bottom:top fov][x, y, z]
# and is turned into a ndarray at the end of calculations.
self.instrument_los = []
# ============= Populating Properties =============
self.__load_data__(source)
# Find for what times should do calculations
if time_range is not None:
if len(time_range) == 0:
self.time_range_bool = np.ones(len(self.TimeStamp), dtype=bool)
elif len(time_range) == 1:
self.time_range_bool = np.zeros(len(self.TimeStamp), dtype=bool)
self.time_range_bool[(np.abs(np.array(self.Time_Hrs) - time_range[0])).argmin()] = 1
elif len(time_range) == 2:
time_range.sort()
self.time_range_bool = np.zeros(len(self.TimeStamp), dtype=bool)
self.time_range_bool[(np.abs(np.array(self.Time_Hrs) - time_range[0])).argmin():
(np.abs(np.array(self.Time_Hrs) - time_range[1])).argmin()] = 1
self.time_range_index = np.where(self.time_range_bool)[0]
self.__get_aircraft_local_vectors__()
self.__get_aircraft_attitude_vectors__()
self.__get_instr_los__()
self.__get_mjds__()
self.__convert_lists_to_ndarrays__()
# This will populate an array of boolean values that filter out the default "start up" points which contain
# no useful information about the flight.
if time_range is not None:
self.good = np.abs(self.aircraft_position_geocentric_xyz[:, 0] - 6378928) > 1000
self.good = np.logical_and(self.good, self.aircraft_position_geocentric_xyz[:, 1] != 0)
self.good = np.logical_and(self.good, self.aircraft_position_geocentric_xyz[:, 2] != 0)
return
# ==================================================================================================================
# -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_Calculation Tools-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
# ==================================================================================================================
def __load_data__(self, source: str) -> None:
"""
This function does the actual loading in, saving, and class property population of the raw
flight data from the url.
:param source: (string) Source of the flight data passed by the constructor.
Original Author: Jeff Langille
Current Author: Daniel Letros
Date Revised: 2017-02-24
"""
real_time = False
flight_retrieve_file = []
if source.lower().find("sim".lower()) > -1:
sim_file_ID_index = source.lower().find("sim".lower()) + 3
flight_retrieve_file = "Flight_Simulation_" + str(source[sim_file_ID_index::]) + ".csv"
if self.__verbose__: print(
"Loading Simulated ER-2 Flight data from " + flight_retrieve_file + "..")
elif source.lower().find("asp-archive.arc.nasa.gov".lower()) > -1:
flight_retrieve_file = "ER2_flight_data.csv"
self.__create_ER2_csv_file__(source, flight_retrieve_file)
else:
# Assuming that it is a one line real time input
real_time = True
if not real_time:
if self.__verbose__: print("Parsing downloaded data..")
self.TimeStamp = self.__flight_data__(flight_retrieve_file, 1) # B
self.Latitude = self.__flight_data__(flight_retrieve_file, 2) # C
self.Longitude = self.__flight_data__(flight_retrieve_file, 3) # D
self.GPS_Altitude_MSL = self.__flight_data__(flight_retrieve_file, 4) # [m] E
self.GPS_Altitude = self.__flight_data__(flight_retrieve_file, 5) # [m] F
self.Pressure_Altitude = self.__flight_data__(flight_retrieve_file, 6) # G
self.RADAR_Altitude = self.__flight_data__(flight_retrieve_file, 7) # [ft] H
self.Ground_Speed = self.__flight_data__(flight_retrieve_file, 8) # [m/s] I
self.True_Air_Speed = self.__flight_data__(flight_retrieve_file, 9) # [m/s] J
self.Indicated_Air_Speed = self.__flight_data__(flight_retrieve_file, 10) # [m/s] K
self.Mach_Number = self.__flight_data__(flight_retrieve_file, 11) # L
self.Vertical_Speed = self.__flight_data__(flight_retrieve_file, 12) # [m/s] M
self.True_Heading = self.__flight_data__(flight_retrieve_file, 13) # [+- 180 deg] N
self.Track_Angle = self.__flight_data__(flight_retrieve_file, 14) # O
self.Drift_Angle = self.__flight_data__(flight_retrieve_file, 15) # [+- 180 deg] P
self.Pitch_Angle = self.__flight_data__(flight_retrieve_file, 16) # [+- 180 deg] Q
self.Roll_Angle = self.__flight_data__(flight_retrieve_file, 17) # [+- 180 deg] R
self.Slip_Angle = self.__flight_data__(flight_retrieve_file, 18) # [deg] S
self.Attack_Angle = self.__flight_data__(flight_retrieve_file, 19) # [deg] T
self.Static_Air_Temp = self.__flight_data__(flight_retrieve_file, 20) # [C] U
self.Dew_Point = self.__flight_data__(flight_retrieve_file, 21) # [C] V
self.Total_Air_Temp = self.__flight_data__(flight_retrieve_file, 22) # [C] W
self.Static_Pressure = self.__flight_data__(flight_retrieve_file, 23) # [mbar] X
self.Dynamic_Pressure = self.__flight_data__(flight_retrieve_file, 24) # [mbar] Y
self.Cabin_Pressure = self.__flight_data__(flight_retrieve_file, 25) # [mbar] Z
self.Wind_Speed = self.__flight_data__(flight_retrieve_file, 26) # [0-132 m/s] AA
self.Wind_Direction = self.__flight_data__(flight_retrieve_file, 27) # AB
self.Vert_Wind_Speed = self.__flight_data__(flight_retrieve_file, 28) # [m/s] AC
self.Solar_Zenith_Angle = self.__flight_data__(flight_retrieve_file, 29) # [deg] AD
self.Aircraft_Sun_Elevation = self.__flight_data__(flight_retrieve_file, 30) # [deg] AE
self.Sun_Azimuth = self.__flight_data__(flight_retrieve_file, 31) # [deg] AF
self.Aircraft_Sun_Azimuth = self.__flight_data__(flight_retrieve_file, 32) # [deg] AG
else:
source_split = np.array(source.split(','))
source_split[source_split == ''] = 0.0
self.TimeStamp = np.array([str(source_split[1])])
self.Latitude = np.array([float(source_split[2])])
self.Longitude = np.array([float(source_split[3])])
self.GPS_Altitude_MSL = np.array([float(source_split[4])]) # [m]
self.GPS_Altitude = np.array([float(source_split[5])]) # [m]
self.Pressure_Altitude = np.array([float(source_split[6])])
self.RADAR_Altitude = np.array([float(source_split[7])]) # [ft]
self.Ground_Speed = np.array([float(source_split[8])]) # [m/s]
self.True_Air_Speed = np.array([float(source_split[9])]) # [m/s]
self.Indicated_Air_Speed = np.array([float(source_split[10])]) # [m/s]
self.Mach_Number = np.array([float(source_split[11])])
self.Vertical_Speed = np.array([float(source_split[12])]) # [m/s]
self.True_Heading = np.array([float(source_split[13])]) # [+- 180 deg]
self.Track_Angle = np.array([float(source_split[14])])
self.Drift_Angle = np.array([float(source_split[15])]) # [+- 180 deg]
self.Pitch_Angle = np.array([float(source_split[16])]) # [+- 180 deg]
self.Roll_Angle = np.array([float(source_split[17])]) # [+- 180 deg]
self.Slip_Angle = np.array([float(source_split[18])]) # [deg]
self.Attack_Angle = np.array([float(source_split[19])]) # [deg]
self.Static_Air_Temp = np.array([float(source_split[20])]) # [C]
self.Dew_Point = np.array([float(source_split[21])]) # [C]
self.Total_Air_Temp = np.array([float(source_split[22])]) # [C]
self.Static_Pressure = np.array([float(source_split[23])]) # [mbar]
self.Dynamic_Pressure = np.array([float(source_split[24])]) # [mbar]
self.Cabin_Pressure = np.array([float(source_split[25])]) # [mbar]
self.Wind_Speed = np.array([float(source_split[26])]) # [0-132 m/s]
self.Wind_Direction = np.array([float(source_split[27])])
self.Vert_Wind_Speed = np.array([float(source_split[28])]) # [m/s]
self.Solar_Zenith_Angle = np.array([float(source_split[29])]) # [deg]
self.Aircraft_Sun_Elevation = np.array([float(source_split[30])]) # [deg]
self.Sun_Azimuth = np.array([float(source_split[31])]) # [deg]
self.Aircraft_Sun_Azimuth = np.array([float(source_split[32])]) # [deg]
self.Time_Hrs = np.zeros(np.array([self.TimeStamp]).size)
for idx in range(len(self.TimeStamp)):
self.Time_Hrs[idx] = float(self.__time_to_hours__(self.TimeStamp[idx]))
return
def __time_to_hours__(self, time_stamp: str) -> float:
"""
This function will convert the ER-2 date time stamp to hours.
:param time_stamp: (string) The ER-2 time stamp.
:return: (float) The converted time stamp.
Original Author: Jeff Langille
Current Author: Daniel Letros
Date Revised: 2017-02-24
"""
try:
dateyy, time = time_stamp.split('T')
hour_min_sec, milli = time.split('.')
hours, minutes, seconds = hour_min_sec.split(':')
return int(hours) + int(minutes) / 60 + int(seconds) / 3600.0
except:
return -99
def __flight_data__(self, filename: str, column_num: int) -> list:
"""
This function reads a column of the csv data file containing the ER-2 Flight Data.
:param filename: (string) The file name of the csv data file.
:param column_num: (int) The column number of the specific data to be retrieved.
:return: (list) The flight data from the column.
Original Author: Jeff Langille
Current Author: Daniel Letros
Date Revised: 2017-02-24
"""
with open(filename) as csvfile:
readCSV = csv.reader(csvfile, delimiter=',')
data = []
for column in readCSV:
values = column[column_num]
data.append(values)
data_out = np.zeros(np.array([data]).size)
if column_num > 1:
for idx in range(np.array([data]).size):
try:
data_out[idx] = float(data[idx])
except:
data_out[idx] = 0
else:
for idx in range(np.array([data]).size):
data_out = data
return data_out
def __create_ER2_csv_file__(self, url: str, flight_retrieve_file: str) -> None:
"""
This function makes a CSV file and writes data from the URL to the CSV file for later parsing.
:param url: (string) The url link to a specific flight.
:param flight_retrieve_file: (string) Name of the CSV file to save the data in.
Original Author: Jeff Langille
Current Author: Daniel Letros
Date Revised: 2017-02-24
"""
if self.__verbose__: print("Downloading ER-2 Flight data from " + str(url) + " to " + str(flight_retrieve_file) + "..")
f = urls.urlopen(url, context=ssl._create_unverified_context()) # To get rid of unverified error, probably not great.
data = f.read()
with open(flight_retrieve_file, "wb") as code:
code.write(data)
return
def __get_aircraft_local_vectors__(self):
"""
This function populates self.aircraft_position_geocentric_xyz, self.aircraft_local_up_geocentric_xyz,
self.aircraft_local_down_geocentric_xyz, self.aircraft_local_north_geocentric_xyz,
self.aircraft_local_south_geocentric_xyz, self.aircraft_local_east_geocentric_xyz,
and self.aircraft_local_west_geocentric_xyz.
Author: Daniel Letros
Date Revised: 2017-03-01
"""
self.__get_aircraft_position_geocentric_xyz__()
self.__get_aircraft_local_up_down_geocentric_xyz__()
self.__get_aircraft_local_north_south_geocentric_xyz__()
self.__get_aircraft_local_east_west_geocentric_xyz__()
def __get_aircraft_attitude_vectors__(self):
"""
This function populates self.aircraft_attitude_up_geocentric_xyz, self.aircraft_attitude_down_geocentric_xyz,
self.aircraft_attitude_forward_geocentric_xyz, self.aircraft_attitude_backward_geocentric_xyz,
self.aircraft_attitude_left_geocentric_xyz, and self.aircraft_attitude_right_geocentric_xyz.
Author: Daniel Letros
Date Revised: 2017-03-01
"""
if self.__verbose__: print("Finding aircraft attitude orientation..")
for time_idx in range(len(self.time_range_index)):
# For start assume the aircraft is level facing north. Note that all of these get updated in
# self.__do_attitude_rotation__ for each rotation done.
self.aircraft_attitude_forward_geocentric_xyz.append(self.aircraft_local_north_geocentric_xyz[time_idx])
self.aircraft_attitude_backward_geocentric_xyz.append(self.aircraft_local_south_geocentric_xyz[time_idx])
self.aircraft_attitude_left_geocentric_xyz.append(self.aircraft_local_west_geocentric_xyz[time_idx])
self.aircraft_attitude_right_geocentric_xyz.append(self.aircraft_local_east_geocentric_xyz[time_idx])
self.aircraft_attitude_up_geocentric_xyz.append(self.aircraft_local_up_geocentric_xyz[time_idx])
self.aircraft_attitude_down_geocentric_xyz.append(self.aircraft_local_down_geocentric_xyz[time_idx])
# Rotate aircraft to face true heading. To match with definition of units do so by rotating about up
# vector by the true heading.
self.__do_attitude_rotation__(self.aircraft_attitude_up_geocentric_xyz[time_idx],
self.True_Heading[self.time_range_index[time_idx]]*np.pi/180,
time_idx)
# Rotate aircraft for pitch. To match with definition of units do so by rotating about the left vector
# by the pitch angle.
self.__do_attitude_rotation__(self.aircraft_attitude_left_geocentric_xyz[time_idx],
self.Pitch_Angle[self.time_range_index[time_idx]]*np.pi/180,
time_idx)
# Rotate aircraft for roll. To match with definition of units do so by rotating about the backward vector
# by the roll angle.
self.__do_attitude_rotation__(self.aircraft_attitude_backward_geocentric_xyz[time_idx],
self.Roll_Angle[self.time_range_index[time_idx]]*np.pi/180,
time_idx)
return
def __get_aircraft_position_geocentric_xyz__(self) -> None:
"""
This function gets the aircraft position in geocentric x y and z coordinates from the raw flight data.
Inputs and outputs are all class internal. Uses self.Latitude, self.Longitude, and self.GPS_Altitude_MSL to
populate self.aircraft_position_geocentric_xyz.
Author: Daniel Letros
Date Revised: 2017-02-24
"""
if self.__verbose__: print("Finding aircraft geocentric coordinates..")
latitude_rad = self.Latitude * (np.pi/180)
longitude_rad = self.Longitude * (np.pi/180)
altitude = self.GPS_Altitude_MSL
num_index = len(self.time_range_index)
# # ============================================================================================================
# # Can use this code if willing to include sasktranif package. For WGS84 this is in near
# # perfect (down to about 10^-8) agreement with the code below.
# geoid = saskif.ISKGeodetic()
# geoid.SetProperty("WGS84", 1)
# for idx in range(num_index):
# geoid.SetLocationLatLonAlt(latitude_rad[idx] * 180 / np.pi, longitude_rad[idx] * 180 / np.pi, altitude[idx])
# xyz_list = geoid.GetLocationXYZ()
# self.aircraft_position_geocentric_xyz.append(xyz_list)
# # ============================================================================================================
for idx in range(num_index):
self.aircraft_position_geocentric_xyz.append(
self.__convert_lat_lon_alt_to_geocentric_xyz__(latitude_rad[self.time_range_index[idx]],
longitude_rad[self.time_range_index[idx]],
altitude[self.time_range_index[idx]]))
return
def __get_aircraft_local_up_down_geocentric_xyz__(self):
"""
This function gets the normalized aircraft local up and down vectors in geocentric x y and z coordinates
from the raw flight data.
Inputs and outputs are all class internal. Uses self.Latitude, self.Longitude, and self.GPS_Altitude_MSL to
populate self.aircraft_local_up_geocentric_xyz and self.aircraft_local_down_geocentric_xyz.
Author: Daniel Letros
Date Revised: 2017-02-24
"""
if self.__verbose__: print("Finding aircraft geocentric local up/down..")
latitude_rad = self.Latitude * (np.pi / 180)
longitude_rad = self.Longitude * (np.pi / 180)
altitude = self.GPS_Altitude_MSL
altitude_higher = self.GPS_Altitude_MSL + 10000
num_index = len(self.time_range_index)
# # ============================================================================================================
# # Can use this code if willing to include sasktranif package. For WGS84 this is in near
# # perfect (down to about 10^-13) agreement with the code below.
# geoid = saskif.ISKGeodetic()
# geoid.SetProperty("WGS84", 1)
# for idx in range(num_index):
# geoid.SetLocationLatLonAlt(latitude_rad[idx] * 180 / np.pi, longitude_rad[idx] * 180 / np.pi, altitude[idx])
# up_list = geoid.GetLocalUp()
# down_list = [-1*up_list[0], -1*up_list[1], -1*up_list[2]]
# self.aircraft_local_up_geocentric_xyz.append(up_list)
# self.aircraft_local_down_geocentric_xyz.append(down_list)
# # ============================================================================================================
for idx in range(num_index):
v_2 = self.__convert_lat_lon_alt_to_geocentric_xyz__(latitude_rad[self.time_range_index[idx]],
longitude_rad[self.time_range_index[idx]],
altitude_higher[self.time_range_index[idx]])
v_1 = self.__convert_lat_lon_alt_to_geocentric_xyz__(latitude_rad[self.time_range_index[idx]],
longitude_rad[self.time_range_index[idx]],
altitude[self.time_range_index[idx]])
up = self.__get_normalized_difference_vector__(v_2, v_1)
down = [-1*up[0], -1*up[1], -1*up[2]]
self.aircraft_local_up_geocentric_xyz.append(up)
self.aircraft_local_down_geocentric_xyz.append(down)
return
def __get_aircraft_local_north_south_geocentric_xyz__(self):
"""
This function gets the normalized aircraft local north and south vector in geocentric x y and z coordinates
from the raw flight data.
Inputs and outputs are all class internal. Uses self.Latitude, self.Longitude, and self.GPS_Altitude_MSL to
populate self.aircraft_local_north_geocentric_xyz and self.aircraft_local_south_geocentric_xyz.
Author: Daniel Letros
Date Revised: 2017-02-24
"""
if self.__verbose__: print("Finding aircraft geocentric local north/south..")
latitude_rad = self.Latitude * (np.pi / 180)
latitude_rad_norther = (self.Latitude + 0.000001) * (np.pi / 180)
longitude_rad = self.Longitude * (np.pi / 180)
altitude = self.GPS_Altitude_MSL
num_index = len(self.time_range_index)
# # ============================================================================================================
# # Can use this code if willing to include sasktranif package. For WGS84 this is in near
# # perfect (down to about 10^-8) agreement with the code below.
# geoid = saskif.ISKGeodetic()
# geoid.SetProperty("WGS84", 1)
# for idx in range(num_index):
# geoid.SetLocationLatLonAlt(latitude_rad[idx] * 180 / np.pi, longitude_rad[idx] * 180 / np.pi, altitude[idx])
# south_list = geoid.GetLocalSouth()
# north_list = [-1*south_list[0], -1*south_list[1], -1*south_list[2]]
# self.aircraft_local_south_geocentric_xyz.append(south_list)
# self.aircraft_local_north_geocentric_xyz.append(north_list)
# # ============================================================================================================
for idx in range(num_index):
v_2 = self.__convert_lat_lon_alt_to_geocentric_xyz__(latitude_rad_norther[self.time_range_index[idx]],
longitude_rad[self.time_range_index[idx]],
altitude[self.time_range_index[idx]])
v_1 = self.__convert_lat_lon_alt_to_geocentric_xyz__(latitude_rad[self.time_range_index[idx]],
longitude_rad[self.time_range_index[idx]],
altitude[self.time_range_index[idx]])
north = self.__get_normalized_difference_vector__(v_2, v_1)
south = [-1 * north[0], -1 * north[1], -1 * north[2]]
self.aircraft_local_north_geocentric_xyz.append(north)
self.aircraft_local_south_geocentric_xyz.append(south)
return
def __get_aircraft_local_east_west_geocentric_xyz__(self):
"""
This function gets the normalized aircraft local east and west vector in geocentric x y and z coordinates
from the raw flight data.
Inputs and outputs are all class internal. Uses self.Latitude, self.Longitude, and self.GPS_Altitude_MSL to
populate self.aircraft_local_east_geocentric_xyz and self.aircraft_local_west_geocentric_xyz.
Author: Daniel Letros
Date Revised: 2017-02-24
"""
if self.__verbose__: print("Finding aircraft geocentric local east/west..")
latitude_rad = self.Latitude * (np.pi / 180)
longitude_rad = self.Longitude * (np.pi / 180)
longitude_rad_easter = (self.Longitude + 0.000001) * (np.pi / 180)
altitude = self.GPS_Altitude_MSL
num_index = len(self.time_range_index)
# # ============================================================================================================
# # Can use this code if willing to include sasktranif package. For WGS84 this is in near
# # perfect (down to about 10^-8) agreement with the code below.
# geoid = saskif.ISKGeodetic()
# geoid.SetProperty("WGS84", 1)
# for idx in range(num_index):
# geoid.SetLocationLatLonAlt(latitude_rad[idx] * 180 / np.pi, longitude_rad[idx] * 180 / np.pi, altitude[idx])
# west_list = geoid.GetLocalWest()
# east_list = [-1*west_list[0], -1*west_list[1], -1*west_list[2]]
# self.aircraft_local_west_geocentric_xyz.append(west_list)
# self.aircraft_local_east_geocentric_xyz.append(east_list)
# # ============================================================================================================
for idx in range(num_index):
v_2 = self.__convert_lat_lon_alt_to_geocentric_xyz__(latitude_rad[self.time_range_index[idx]],
longitude_rad_easter[self.time_range_index[idx]],
altitude[self.time_range_index[idx]])
v_1 = self.__convert_lat_lon_alt_to_geocentric_xyz__(latitude_rad[self.time_range_index[idx]],
longitude_rad[self.time_range_index[idx]],
altitude[self.time_range_index[idx]])
east = self.__get_normalized_difference_vector__(v_2, v_1)
west = [-1 * east[0], -1 * east[1], -1 * east[2]]
self.aircraft_local_east_geocentric_xyz.append(east)
self.aircraft_local_west_geocentric_xyz.append(west)
return
def __convert_lat_lon_alt_to_geocentric_xyz__(self, latitude: float, longitude: float, altitude: float) -> list:
"""
This function turns geodetic latitude, longitude and altitude into geocentric x y and z coordinates.
This function uses the WGS84 ellipsoid and makes no geoid correction.
:param latitude: (float) Geodetic latitude in radians.
:param longitude: (float) Geodetic longitude in radians.
:param altitude: (float) Altitude above mean sea level in meters.
:return: (list) The geocentric [x, y, z] coordinates
Author: Daniel Letros
Date Revised: 2017-02-24
"""
semi_major = 6378137 # Semi-major [m]
flattening = 1 / 298.257223563 # Flattening
c = 1 / np.sqrt(np.cos(latitude) ** 2 + (1 - flattening) ** 2 * np.sin(latitude) ** 2)
s = (1 - flattening) ** 2 * c
x = (semi_major * c + altitude) * np.cos(latitude) * np.cos(longitude)
y = (semi_major * c + altitude) * np.cos(latitude) * np.sin(longitude)
z = (semi_major * s + altitude) * np.sin(latitude)
return [x, y, z]
def __get_normalized_difference_vector__(self, v_2: list, v_1: list) -> list:
"""
This function calculates the normalized vector which points from v_1 towards v_2 and returns it.
:param v_2: (list) A vector in the form of [x, y, z]
:param v_1: (list) A vector in the form of [x, y, z]
:return: (list) A normalized vector in the form of [x, y, z]
Author: Daniel Letros
Date Revised: 2017-02-24
"""
vector_x = v_2[0] - v_1[0]
vector_y = v_2[1] - v_1[1]
vector_z = v_2[2] - v_1[2]
magnitude = np.sqrt(vector_x**2 + vector_y**2 + vector_z**2)
return [vector_x/magnitude, vector_y/magnitude, vector_z/magnitude]
def __get_rotation_about_arbitrary_axis__(self, rotation_axis: list, theta: float) -> np.matrix:
"""
This function makes a 3d rotation matrix which will rotate a vector by an angle about a specified axis.
:param rotation_axis: (list) The unit axis the rotation matrix should rotate about [x, y, z]
:param theta: (float) The angle in radians the rotation matrix should do
:return: (np.matrix) rotation matrix in the form of [[r11, r12, r13], [r21, r22, r23], [r31, r32, r33]]
Author: Daniel Letros
Date Revised: 2017-03-01
"""
u = np.array([rotation_axis[0], rotation_axis[1], rotation_axis[2]])
r11 = np.cos(theta) + (u[0] ** 2) * (1 - np.cos(theta))
r12 = u[0] * u[1] * (1 - np.cos(theta)) - u[2] * np.sin(theta)
r13 = u[0] * u[2] * (1 - np.cos(theta)) + u[1] * np.sin(theta)
r21 = u[1] * u[0] * (1 - np.cos(theta)) + u[2] * np.sin(theta)
r22 = np.cos(theta) + (u[1] ** 2) * (1 - np.cos(theta))
r23 = u[1] * u[2] * (1 - np.cos(theta)) - u[0] * np.sin(theta)
r31 = u[2] * u[0] * (1 - np.cos(theta)) - u[1] * np.sin(theta)
r32 = u[2] * u[1] * (1 - np.cos(theta)) + u[0] * np.sin(theta)
r33 = np.cos(theta) + (u[2] ** 2) * (1 - np.cos(theta))
return np.matrix([[r11, r12, r13], [r21, r22, r23], [r31, r32, r33]])
def __do_attitude_rotation__(self, axis: list, theta: float, time_idx: int):
"""
This function will apply a rotation about the inputted axis by angle theta to all aircraft attitude vectors.
:param axis: (list) The unit axis the rotation matrix should rotate about [x, y, z]
:param theta: (float) The angle in radians the rotation matrix should do
:param time_idx: (int) The time index to do the rotation at
Author: Daniel Letros
Date Revised: 2017-03-01
"""
forward = self.aircraft_attitude_forward_geocentric_xyz[time_idx]
backward = self.aircraft_attitude_backward_geocentric_xyz[time_idx]
left = self.aircraft_attitude_left_geocentric_xyz[time_idx]
right = self.aircraft_attitude_right_geocentric_xyz[time_idx]
up = self.aircraft_attitude_up_geocentric_xyz[time_idx]
down = self.aircraft_attitude_down_geocentric_xyz[time_idx]
rotation_matrix = self.__get_rotation_about_arbitrary_axis__(axis, theta)
forward *= rotation_matrix
backward *= rotation_matrix
left *= rotation_matrix
right *= rotation_matrix
up *= rotation_matrix
down *= rotation_matrix
# Turn back into lists
self.aircraft_attitude_forward_geocentric_xyz[time_idx] = [forward.item(0),
forward.item(1),
forward.item(2)]
self.aircraft_attitude_backward_geocentric_xyz[time_idx] = [backward.item(0),
backward.item(1),
backward.item(2)]
self.aircraft_attitude_left_geocentric_xyz[time_idx] = [left.item(0),
left.item(1),
left.item(2)]
self.aircraft_attitude_right_geocentric_xyz[time_idx] = [right.item(0),
right.item(1),
right.item(2)]
self.aircraft_attitude_up_geocentric_xyz[time_idx] = [up.item(0),
up.item(1),
up.item(2)]
self.aircraft_attitude_down_geocentric_xyz[time_idx] = [down.item(0),
down.item(1),
down.item(2)]
return
def __get_instr_los__(self) -> None:
"""
This function will use the bore sight information given in the constructor of the class and produce line of
sight vectors for the instrument. This is done by rotating the forward aircraft vector to produce a line of
sight at the top of the field of view and at the bottom of the field of view with other lines equally
distributed in between to produce a total of self.num_los lines of slight.
Author: Daniel Letros
Date Revised: 2017-03-31
"""
if self.__verbose__: print("Finding instrument lines of sight..")
# Last LOS is equal to forward vector, First is equal to fov angle lower then forward vector (not including
# bore offsets).
fov_angles = np.flipud(np.linspace(0, self.vert_fov * np.pi/180, self.num_los)) # flipped so the "bottom" LOS is
# first in order. Done so
# tangent altitudes are ascending
for idx_time in range(len(self.time_range_index)): # For each point in time..
# Fan the forward vector down to get the LOS
self.instrument_los.append([])
for idx_fov_angle in range(len(fov_angles)):
# Account for vertical bore offset when finding the rotation matrix for calculating the current los
rot_matrix_vert = self.__get_rotation_about_arbitrary_axis__(
self.aircraft_attitude_right_geocentric_xyz[idx_time],
fov_angles[idx_fov_angle] - ((self.bore_offset_top + self.bore_offset_vert) * np.pi/180))
# Account for horizontal bore offset in current los
rot_matrix_horz = self.__get_rotation_about_arbitrary_axis__(
self.aircraft_attitude_up_geocentric_xyz[idx_time],
self.bore_offset_horz * np.pi / 180)
los = self.aircraft_attitude_forward_geocentric_xyz[idx_time] * rot_matrix_vert * rot_matrix_horz
self.instrument_los[idx_time].append([los.item(0), los.item(1), los.item(2)])
return
def __get_mjds__(self) -> None:
"""
This function converts all time stamps into mjds.
Author: Daniel Letros
Date Revised: 2017-04-19
"""
time = ast_time(np.array(self.TimeStamp)[self.time_range_bool], format='isot', scale='utc', precision=9)
self.mjd = time.mjd
def __convert_lists_to_ndarrays__(self) -> None:
"""
This function will take all the calculated data describing aircraft position and orientation and transform them
from a list of lists to a numpy array.
Author: Daniel Letros
Date Revised: 2017-03-24
"""
self.aircraft_position_geocentric_xyz = np.array(self.aircraft_position_geocentric_xyz)
self.aircraft_local_up_geocentric_xyz = np.array(self.aircraft_local_up_geocentric_xyz)
self.aircraft_local_down_geocentric_xyz = np.array(self.aircraft_local_down_geocentric_xyz)
self.aircraft_local_north_geocentric_xyz = np.array(self.aircraft_local_north_geocentric_xyz)
self.aircraft_local_south_geocentric_xyz = np.array(self.aircraft_local_south_geocentric_xyz)
self.aircraft_local_east_geocentric_xyz = np.array(self.aircraft_local_east_geocentric_xyz)
self.aircraft_local_west_geocentric_xyz = np.array(self.aircraft_local_west_geocentric_xyz)
self.aircraft_attitude_up_geocentric_xyz = np.array(self.aircraft_attitude_up_geocentric_xyz)
self.aircraft_attitude_down_geocentric_xyz = np.array(self.aircraft_attitude_down_geocentric_xyz)
self.aircraft_attitude_forward_geocentric_xyz = np.array(self.aircraft_attitude_forward_geocentric_xyz)
self.aircraft_attitude_backward_geocentric_xyz = np.array(self.aircraft_attitude_backward_geocentric_xyz)
self.aircraft_attitude_left_geocentric_xyz = np.array(self.aircraft_attitude_left_geocentric_xyz)
self.aircraft_attitude_right_geocentric_xyz = np.array(self.aircraft_attitude_right_geocentric_xyz)
self.instrument_los = np.array(self.instrument_los)
# ==================================================================================================================
# -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_Visualization Tools-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
# ==================================================================================================================
[docs] def plot_flight_path(self, fig_num: int=1, plot_wgs84: bool=False, flat=False) -> None:
"""
This function makes a plot of the aircraft flight path. If the plot is to be made on the WGS84 ellipsoid then
all coordinates will be in the geocentric system.
If the plot is NOT to be made on the ellipsoid then the path of the aircraft can be plotted in two
different ways, flattened and un-flattened.
The un-flattened method plots the aircraft as tho it took off from latitude = 0 and longitude = 0.
Position is then found under these conditions in the geocentric system which then undergoes a coordinate
rotation such that aircraft up is along the z axis. The z axis is then scaled such that sea level = 0 meters.
Note that due to the varying nature of the curvature of the WGS84 ellipsoid that the numerical results will be
inaccurate compared to the true geocentric of the aircraft position, however, this provides an excellent
visualization of the aircraft flight path with respect to the curvature of the Earth and time.
The flattened method works just like the un-flattened except that the z-axis is replaced with the mean sea
level altitude provided by the flight GPS.
:param fig_num: (int, optional, default=1) figure number to make plot on
:param plot_wgs84: (bool, optional, default=False) If true will also plot the WGS84 ellipsoid so you can
visualize the path with respect to the globe. Flat input is ignored, plot is in true
geocentric coordinates.
:param flat: (bool, optional, default=False) If true will plot the flight path without any curvature of the
Earth.
Author: Daniel Letros
Date Revised: 2017-03-24
"""
if self.__verbose__: print("Plotting flight path..")
fig = mplot.figure(num=fig_num)
ax = fig.add_subplot(111, projection='3d')
if plot_wgs84:
self.__plot_wsg84_elip__()
ax.plot(self.aircraft_position_geocentric_xyz[self.good, 0],
self.aircraft_position_geocentric_xyz[self.good, 1],
self.aircraft_position_geocentric_xyz[self.good, 2])
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_zlabel('z [m]')
else:
# EAFB is lat, lon, alt of 34.923524, -117.909190, 0. Treat this as the zero point. Find lat, lon, alt
# change that the aircraft makes with respect to this point. Now have the lat, lon, and alt of the aircraft
# as tho if took off from a lat, lon, alt of 0, 0, 0.
lat_normed = (self.Latitude[self.time_range_bool] * np.pi / 180) - (34.923524 * np.pi / 180)
lon_normed = (self.Longitude[self.time_range_bool] * np.pi / 180) - (-117.909190 * np.pi / 180)
alt_normed = self.GPS_Altitude_MSL[self.time_range_bool]
geocentric_coordinates_respect_to_norm = self.__convert_lat_lon_alt_to_geocentric_xyz__(lat_normed,
lon_normed,
alt_normed)
geocentric_coordinates_respect_to_norm = np.array(geocentric_coordinates_respect_to_norm)
# Want to plot our x axis as our z axis such that "up" is along z.
if flat:
self.__plot_line_colored_in_time__(-geocentric_coordinates_respect_to_norm[2, :],
geocentric_coordinates_respect_to_norm[1, :],
self.GPS_Altitude_MSL[self.time_range_bool],
self.Time_Hrs[self.time_range_bool])
else:
geocentric_coordinates_level = self.__convert_lat_lon_alt_to_geocentric_xyz__(0, 0, 0)
self.__plot_line_colored_in_time__(-geocentric_coordinates_respect_to_norm[2, :],
geocentric_coordinates_respect_to_norm[1, :],
geocentric_coordinates_respect_to_norm[0, :] -
geocentric_coordinates_level[0],
self.Time_Hrs[self.time_range_bool])
return
[docs] def plot_instantaneous_attitude_solution(self, time_index: int, plot_wgs84: bool=False, fig_num: int=1,
scale_factor: int=1000000) -> None:
"""
This function makes a plot of the attitude of the aircraft in geocentric coordinates for a given point in time.
:param time_index: (int) Index of the time you wish to plot
:param plot_wgs84: (bool, optional, default=False) If true will also plot the WGS84 ellipsoid so you can
visualize the attitude with respect to the globe.
:param fig_num: (int, optional, default=1) figure number to make plot on
:param scale_factor: (int, optional, default=1000000) Used to scale the vectors to be visible on the
WGS84 ellipsoid.
Author: Daniel Letros
Date Revised: 2017-03-01
"""
if self.__verbose__: print("Plotting instantaneous attitude..")
fig = mplot.figure(num=fig_num)
ax = fig.add_subplot(111, projection='3d')
if plot_wgs84:
self.__plot_wsg84_elip__(ax)
# Plot aircraft position and local vectors
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0]],
[self.aircraft_position_geocentric_xyz[time_index][1]],
[self.aircraft_position_geocentric_xyz[time_index][2]],
markerfacecolor='black', markeredgecolor='black', marker='o', markersize=10, alpha=0.6)
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.aircraft_local_up_geocentric_xyz[time_index][0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.aircraft_local_up_geocentric_xyz[time_index][1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.aircraft_local_up_geocentric_xyz[time_index][2] * scale_factor],
label="Local Up", color='darkgoldenrod')
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.aircraft_local_north_geocentric_xyz[time_index][0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.aircraft_local_north_geocentric_xyz[time_index][1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.aircraft_local_north_geocentric_xyz[time_index][2] * scale_factor],
label="Local North", color='firebrick')
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.aircraft_local_west_geocentric_xyz[time_index][0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.aircraft_local_west_geocentric_xyz[time_index][1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.aircraft_local_west_geocentric_xyz[time_index][2] * scale_factor],
label="Local West", color='green')
# Plot aircraft attitude vectors
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.aircraft_attitude_up_geocentric_xyz[time_index][0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.aircraft_attitude_up_geocentric_xyz[time_index][1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.aircraft_attitude_up_geocentric_xyz[time_index][2] * scale_factor],
label="Up", color='yellow')
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.aircraft_attitude_forward_geocentric_xyz[time_index][0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.aircraft_attitude_forward_geocentric_xyz[time_index][1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.aircraft_attitude_forward_geocentric_xyz[time_index][2] * scale_factor],
label="Forward", color='red')
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.aircraft_attitude_left_geocentric_xyz[time_index][0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.aircraft_attitude_left_geocentric_xyz[time_index][1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.aircraft_attitude_left_geocentric_xyz[time_index][2] * scale_factor],
label="Left", color='lawngreen')
# Plot instrument lines of sight
for idx in range(self.num_los):
if idx == 0:
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.instrument_los[time_index][idx, 0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.instrument_los[time_index][idx, 1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.instrument_los[time_index][idx, 2] * scale_factor],
label="los", color='black', linestyle='--', alpha=0.5)
else: # don't label it
mplot.plot([self.aircraft_position_geocentric_xyz[time_index][0],
self.aircraft_position_geocentric_xyz[time_index][0] +
self.instrument_los[time_index][idx, 0] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][1],
self.aircraft_position_geocentric_xyz[time_index][1] +
self.instrument_los[time_index][idx, 1] * scale_factor],
[self.aircraft_position_geocentric_xyz[time_index][2],
self.aircraft_position_geocentric_xyz[time_index][2] +
self.instrument_los[time_index][idx, 2] * scale_factor],
color='black', linestyle='--', alpha=0.5)
mplot.title(str(self.TimeStamp[time_index]) +
" True Heading: " + str(self.True_Heading[self.time_range_index[time_index]]) + " [deg],"
" Pitch Angle: " + str(self.Pitch_Angle[self.time_range_index[time_index]]) + " [deg],"
" Roll Angle: " + str(self.Roll_Angle[self.time_range_index[time_index]]) + " [deg]")
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
if not plot_wgs84:
ax.set_aspect('equal')
mplot.legend()
return
def __plot_line_colored_in_time__(self, x: np.ndarray, y: np.ndarray, z: np.ndarray, t: np.ndarray) -> None:
"""
This function will plot the x, y, z data on the current figure with different colors based on their indices.
The different color segments will then be labeled based on the data values in t.
:param x: (1d ndarray of floats) data to plot on x axis
:param y: (1d ndarray of floats) data to plot on y axis
:param z: (1d ndarray of floats) data to plot on z axis
:param t: (1d ndarray of floats) data to label color segments
Author: Daniel Letros
Date Revised: 2017-03-30
"""
color_chunks = 12
# Use jet cmap to get colors with respect to indices in time
cmap = mplot.cm.Paired(np.linspace(0, 1, color_chunks))
points_in_chunk = np.floor(np.floor(len(t))/color_chunks)
ax = mplot.gca()
start_of_chunk_index = 0
end_of_chunk_index = points_in_chunk
for i in range(color_chunks):
if (end_of_chunk_index + points_in_chunk) > len(t):
# going to fall short of coloring all points (remainder is left) so add extra points to last chunk
ax.plot(x[start_of_chunk_index::],
y[start_of_chunk_index::],
z[start_of_chunk_index::],
color=cmap[i],
label=str(np.round(t[start_of_chunk_index], 2)) + " : " + str(
np.round(t[-1], 2)))
else:
ax.plot(x[start_of_chunk_index:end_of_chunk_index],
y[start_of_chunk_index:end_of_chunk_index],
z[start_of_chunk_index:end_of_chunk_index],
color=cmap[i],
label=str(np.round(t[start_of_chunk_index], 2)) + " : " + str(np.round(t[end_of_chunk_index], 2)))
start_of_chunk_index += points_in_chunk
end_of_chunk_index += points_in_chunk
ax.set_xlabel('South (+) [m]')
ax.set_ylabel('East (+) [m]')
ax.set_zlabel('Up (+) [m]')
mplot.legend()
return
def __plot_wsg84_elip__(self) -> None:
"""
This function will plot the WGS84 ellipsoid on the current figure
Author: Daniel Letros
Date Revised: 2017-03-24
"""
# Plot the ellipsoid as well as the position of the university of saskatchewan and edwards air force base.
# usask and EAFB are taken to be at sea level for visualization purposes.
# -------------------------------------Make and plot WGS84 ellipsoid----------------------------------------
semi_major = 6378137 # Semi-major [m]
flattening = 1 / 298.257223563 # Flattening
a = semi_major
b = a * (1 - flattening)
c = a ** 2 / b
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
WGS84_ellipsoid_x = np.outer(np.cos(u), np.sin(v)) / np.sqrt(1 / (a ** 2))
WGS84_ellipsoid_y = np.outer(np.sin(u), np.sin(v)) / np.sqrt(1 / (b ** 2))
WGS84_ellipsoid_z = np.outer(np.ones_like(u), np.cos(v)) / np.sqrt(1 / (c ** 2))
ax = mplot.gca()
ax.plot_surface(WGS84_ellipsoid_x, WGS84_ellipsoid_y, WGS84_ellipsoid_z,
color='b', linewidth=0, rstride=4, cstride=4)
# ----------------------------------------------------------------------------------------------------------
# ------------------------------------Make and plot Usask and EAFB------------------------------------------
usask_list = self.__convert_lat_lon_alt_to_geocentric_xyz__(52.131459 * np.pi / 180, -106.632935 * np.pi / 180,
0)
EAFB_list = self.__convert_lat_lon_alt_to_geocentric_xyz__(34.923524 * np.pi / 180, -117.909190 * np.pi / 180,
0)
mplot.plot([usask_list[0]],
[usask_list[1]],
[usask_list[2]],
markerfacecolor='lawngreen', markeredgecolor='lawngreen', marker='d', markersize=10, alpha=0.6)
mplot.plot([EAFB_list[0]],
[EAFB_list[1]],
[EAFB_list[2]],
markerfacecolor='cyan', markeredgecolor='cyan', marker='d', markersize=10, alpha=0.6)
# ----------------------------------------------------------------------------------------------------------
return
# a = ER2_Flight_Data("https://asp-archive.arc.nasa.gov/N809NA/FY2017/2017-03-09/IWG1.09Mar2017-2343",
# verbose=True, vert_fov=4, bore_offset_vert=0, bore_offset_horz=0, bore_offset_top=2, num_los=10)
# # a.plot_flight_path(fig_num=1, plot_wgs84=False, flat=False)
# a.plot_instantaneous_attitude_solution(12646, scale_factor=1, fig_num=2)
# mplot.show()