use data from a vector to locate a row in a matrix and extract the data from the row.

5 次查看（过去 30 天）
Terry Poole
Terry Poole 2022-2-28
回答： Gyan Vaibhav 2023-11-18
So I have a 902x1 vector of unique time stamps, and I have 9 different 901x6 arrays of data (column 1 is the time stamp for each row of data)
I need to take the time stamp from the unique_time vector and find the timestamp in each 901x6 array that matches it and extract those rows of data.
i've tried for loops and if statements, but none work, I have to admit this level of complexity is above my head (programatically speaking).
The part of the code in question starts at line 70, there are a lot of functions involved in the script as a whole so i'll add them as well.
In my mind it should go something like
for i = 1:unique_time
for j = 1:9
if unique_time(i) == mat(j)(:, 1)
Data_mat = [mat(j)(:, 2); mat(j)(:, 3); mat(j)(:, 4)];
else
Data_mat = [0, 0, 0];
end
end
end
But this fails miserably. I'm sure it's some trick with the ismember function that I'm just not experienced enough to figure out.
Any help/advise/point in a general direction would be greatly appreciated.
clear all; close all; clc;
%% Setup WGS84 parameters
WGS.EQUATORIAL_RADIUS = 6378137;
WGS.FLATTENING = 1/298.257223563;
WGS.EARTH_ROTATION = 7292115 * 10^(-11); % rad/s
WGS.EARTH_GRAVITIONAL_CONSTANT = 3.9860050 * 10^(14); % m^3/s^2
WGS.ECCENTRICITY = sqrt(WGS.FLATTENING * (2-WGS.FLATTENING));
WGS.ECCENTRICITY_SQUARED = WGS.FLATTENING * (2-WGS.FLATTENING);
WGS.POLAR_RADIUS = WGS.EQUATORIAL_RADIUS * (1-WGS.FLATTENING);
%% Compute average base position
fileID = fopen('data_base\ecef_rx0.txt');
base = textscan(fileID, '%f %d %f %f %f %f %f %f %f %f %f %f %f %f %d %d %d %d');
fclose(fileID);
X = mean(base{1, 3});
Y = mean(base{1, 4});
Z = mean(base{1, 5});
xyz = [X, Y, Z];
base_lla = ecef2lla(xyz, WGS);
%% Compute Sattelite Line of Sight Vectors
% Read in icp satelite data
sat = string({'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat2.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat4.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat5.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat9.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat10.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat12.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat17.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat23.txt';
'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat25.txt'});
% Read in broadcast file
ephemeris = read_rinex_nav('brdc2930.11n');
% Assign each satelite a number
sat_num = 1:1:size(sat, 1);
% For Loopy stuff to find Line of Sight Vectors
for i = 1:1:size(sat, 1)
fileID = fopen(sat(i), 'r');
current_sat = sat_num(i);
base_obs = cell2mat(textscan(fileID, '%f %f %f %f %f %f %f %f'));
eph_rel = ephemeris(find(ephemeris(:, 1) == current_sat), :);
los_vec = [];
for j = 1:1:size(base_obs, 1)
[~, eph_rel_line] = min(abs(base_obs(j, 1) - eph_rel(:, 17)));
ecef_pos = gpseph2ecef(base_obs(j, 1), eph_rel(eph_rel_line, :));
sv_obs = [0, 0, 0, 0, 0, 0];
sv_obs(1) = round(base_obs(j, 1));
ned = ecef2ned([ecef_pos(1), ...
ecef_pos(2), ecef_pos(3)], xyz, WGS);
sv_obs(2) = ned(1);
sv_obs(3) = ned(2);
sv_obs(4) = ned(3);
sv_obs(5) = -asin(sv_obs(4)/sqrt(sv_obs(2)^2 + sv_obs(3)^2 ...
+ sv_obs(4)^2))*180/pi;
sv_obs(6) = base_obs(j, 8);
los_vec = vertcat(los_vec, sv_obs);
end
base_data{i} = los_vec;
fclose(fileID);
end
%% Combine 9 different size arrays into 1 single array
%% with zeroes filling the extra values.
mat1 = base_data{1, 1};
mat1(901, 6) = 0;
mat2 = base_data{1, 2};
mat2(901, 6) = 0;
mat3 = base_data{1, 3};
mat3(901, 6) = 0;
mat4 = base_data{1, 4};
mat4(901, 6) = 0;
mat5 = base_data{1, 5};
mat5(901, 6) = 0;
mat6 = base_data{1, 6};
mat6(901, 6) = 0;
mat7 = base_data{1, 7};
mat7(901, 6) = 0;
mat8 = base_data{1, 8};
mat8(901, 6) = 0;
mat9 = base_data{1, 9};
mat9(901, 6) = 0;
%% Extract X, Y, Z coordinate columns into individual arrays
X_mat = [mat1(:, 2), mat2(:, 2), mat3(:, 2), mat4(:, 2), mat5(:, 2)...
mat6(:, 2), mat7(:, 2), mat8(:, 2), mat9(:, 2)];
Y_mat = [mat1(:, 3), mat2(:, 3), mat3(:, 3), mat4(:, 3), mat5(:, 3)...
mat6(:, 3), mat7(:, 3), mat8(:, 3), mat9(:, 3)];
Z_mat = [mat1(:, 4), mat2(:, 4), mat3(:, 4), mat4(:, 4), mat5(:, 4)...
mat6(:, 4), mat7(:, 4), mat8(:, 4), mat9(:, 4)];
%% Find each unique time step
time_mat = [mat1(:, 1), mat2(:, 1), mat3(:, 1), mat4(:, 1), mat5(:, 1)...
mat6(:, 1), mat7(:, 1), mat8(:, 1), mat9(:, 1)];
time_linear = time_mat(:);
time_sort = sort(time_linear);
time_unique = unique(time_sort);
%% Find highest value for X, Y, and Z matrices at each unique time step
Functions used in script:
function lla = ecef2lla(xyz, WGS)
% Converts the WGS Earth-Centered, Earth-Fixed (ECEF) Cartesian
% coordinates to geodetic latitude, longitude, altitude (LLA).
MU = 3.986005 * 10^14; % m^3/s^2
SPEED_OF_LIGHT = 2.99792458 * 10^8; % m/s
EARTH_ROT_RATE = 7.2921151467 * 10^-5; % rad/s
PI = 3.1415926535898;
FLATTENING = 1 / 298.257223563;
ECCENTRICITY = sqrt(FLATTENING * (2 - FLATTENING));
EQ_RAD = 6378137; % m
POL_RAD = EQ_RAD * (1 - FLATTENING);
WGS.EQUATORIAL_RADIUS = 6378137;
WGS.FLATTENING = 1/298.257223563;
WGS.EARTH_ROTATION = 7292115 * 10^(-11); % rad/s
WGS.EARTH_GRAVITIONAL_CONSTANT = 3.9860050 * 10^(14); % m^3/s^2
WGS.ECCENTRICITY = sqrt(WGS.FLATTENING * (2-WGS.FLATTENING));
WGS.ECCENTRICITY_SQUARED = WGS.FLATTENING * (2-WGS.FLATTENING);
WGS.POLAR_RADIUS = WGS.EQUATORIAL_RADIUS * (1-WGS.FLATTENING);
RAD2DEG = 180/pi;
c = 2.99792458 * 10^8; % m/s
x = xyz(1);
y = xyz(2);
z = xyz(3);
% Calculate longitude.
if ((xyz(1) == 0.0) && (xyz(2) == 0.0))
long = 0.0;
else
long = atan2(xyz(2), xyz(1))*RAD2DEG;
end
% Compute altitude and latitude.
if ((xyz(1) == 0.0) && (xyz(2) == 0.0) && ...
(xyz(3) == 0.0))
error('XYZ at center of earth');
else
% Compute altitude.
p = norm([xyz(1) xyz(2)]);
E = sqrt(WGS.EQUATORIAL_RADIUS^2 - WGS.POLAR_RADIUS^2);
F = 54 * (WGS.POLAR_RADIUS*xyz(3))^2;
G = p^2 + (1-WGS.ECCENTRICITY_SQUARED)*xyz(3)^2 ...
- (WGS.ECCENTRICITY*E)^2;
s = (1 + WGS.ECCENTRICITY_SQUARED^2*F*p^2/G^3 ...
+ sqrt(c^2 + 2*c))^(1/3);
P = (F/(3*G^2)) / ((s + (1/s) + 1)^2);
Q = sqrt(1 + 2*WGS.ECCENTRICITY_SQUARED^2*P);
k_1 = -P*WGS.ECCENTRICITY_SQUARED*p / (1 + Q);
k_2 = 0.5*WGS.EQUATORIAL_RADIUS^2*(1 + 1/Q);
k_3 = -(1 - WGS.ECCENTRICITY_SQUARED)*P*z^2 / (Q*(1 + Q));
k_4 = -0.5*P*p^2;
r_0 = k_1 + sqrt(k_2 + k_3 + k_4);
k_5 = (p - WGS.ECCENTRICITY_SQUARED*r_0);
U = sqrt(k_5^2 + xyz(1)^2);
V = sqrt(k_5^2 + (1 - WGS.ECCENTRICITY_SQUARED)*xyz(3)^2);
alt = U*(1 - (WGS.POLAR_RADIUS^2 ...
/(WGS.EQUATORIAL_RADIUS*V)));
% Compute latitude.
z_0 = (WGS.POLAR_RADIUS^2*xyz(3)) ...
/ (WGS.EQUATORIAL_RADIUS*V);
e_p = (WGS.EQUATORIAL_RADIUS/WGS.POLAR_RADIUS) ...
* WGS.ECCENTRICITY;
lat = atan((z + z_0*(e_p)^2)/p) * RAD2DEG;
end
lla = [lat; long; alt];
end
function ephemeris = read_rinex_nav( filename )
% Read the RINEX navigation file. The header is skipped because
% information in it (a0, a1, iono alpha and beta parameters) is not
% currently needed for orbit computation. This can be easily amended to
% include navigation header information by adding lines in the 'while' loop
% where the header is currently skipped.
%
% Input: - filename - enter the filename to be read. If filename
% exists, the orbit will be calculated.
%
% Output: - ephemeris - Output is a matrix with rows for each PRN and
% columns as follows:
%
% col 1: PRN ....... satellite PRN
% col 2: M0 ....... mean anomaly at reference time
% col 3: delta_n ..... mean motion difference
% col 4: e ....... eccentricity
% col 5: sqrt(A) ..... where A is semimajor axis
% col 6: OMEGA ....... LoAN at weekly epoch
% col 7: i0 ....... inclination at reference time
% col 8: omega ....... argument of perigee
% col 9: OMEGA_dot ... rate of right ascension
% col 10: i_dot ....... rate of inclination angle
% col 11: Cuc ....... cosine term, arg. of latitude
% col 12: Cus ....... sine term, arg. of latitude
% col 13: Crc ....... cosine term, radius
% col 14: Crs ....... sine term, radius
% col 15: Cic ....... cosine term, inclination
% col 16: Cis ....... sine term, inclination
% col 17: toe ....... time of ephemeris
% col 18: IODE ....... Issue of Data Ephemeris
% col 19: GPS_wk ....... GPS week
fid = fopen(filename);
if fid == -1
errordlg(['The file ''' filename ''' does not exist.']);
return;
end
% skip through header
end_of_header = 0;
while end_of_header == 0
current_line = fgetl(fid);
if strfind(current_line,'END OF HEADER')
end_of_header=1;
end
end
j = 0;
while feof(fid) ~= 1
j = j+1;
current_line = fgetl(fid);
% parse epoch line (ignores SV clock bias, drift, and drift rate)
line_data = cell2mat(textscan(current_line, '%2.0f %3.0f %3.0f %3.0f %3.0f %3.0f %5.1f %19.12f %19.12f %19.12f'));
Y = line_data(2) + 2000;
M = line_data(3);
D = line_data(4);
H = line_data(5);
min = line_data(6);
sec = line_data(7);
[ gps_week, tow ] = gpsweek(Y, M, D, H, min, sec);
PRN = line_data(1);
af0 = line_data(8);
af1 = line_data(9);
af2 = line_data(10);
% Broadcast orbit line 1
current_line = fgetl(fid);
line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
IODE = line_data(1);
Crs = line_data(2);
delta_n = line_data(3);
M0 = line_data(4);
% Broadcast orbit line 2
current_line = fgetl(fid);
line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
Cuc = line_data(1);
e = line_data(2);
Cus = line_data(3);
sqrtA = line_data(4);
% Broadcast orbit line 3
current_line = fgetl(fid);
line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
toe = line_data(1);
Cic = line_data(2);
OMEGA = line_data (3);
Cis = line_data(4);
% Broadcast orbit line 4
current_line = fgetl(fid);
line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
i0 = line_data(1);
Crc = line_data(2);
omega = line_data (3);
OMEGA_dot = line_data(4);
% Broadcast orbit line 5
current_line = fgetl(fid);
line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
i_dot = line_data(1);
L2_codes = line_data(2);
GPS_wk = line_data (3);
L2_dataflag = line_data(4);
% Broadcast orbit line 6
current_line = fgetl(fid);
line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
SV_acc = line_data(1);
SV_health = line_data(2);
TGD = line_data (3);
IODC = line_data(4);
% Broadcast orbit line 7
current_line = fgetl(fid);
line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
msg_trans_t = line_data(1);
fit_int = line_data(2);
ephemeris(j,:) = [PRN, M0, delta_n, e, sqrtA, OMEGA, i0, omega, OMEGA_dot, i_dot, Cuc, Cus, Crc, Crs, Cic, Cis, toe, IODE, GPS_wk];
end
fclose(fid);
end
function ecef = gpseph2ecef(time, ephemeris)
% Converts a GPS ephemeris to the WGS84 Earth-Centered, Earth-Fixed
% (ECEF) Cartesian coordinates at the specified time.
%
M0 = ephemeris(2);
delta_n = ephemeris(3);
e = ephemeris(4);
sqrtA = ephemeris(5);
OMEGA = ephemeris(6);
i0 = ephemeris(7);
omega = ephemeris(8);
OMEGA_dot = ephemeris(9);
i_dot = ephemeris(10);
Cuc = ephemeris(11);
Cus = ephemeris(12);
Crc = ephemeris(13);
Crs = ephemeris(14);
Cic = ephemeris(15);
Cis = ephemeris(16);
toe = ephemeris(17);
% WGS84 Parameters
MU = 3.986005 * 10^14; % m^3/s^2
EARTH_ROT_RATE = 7.2921151467 * 10^-5; % rad/s
% Compute the semimajor axis.
A = sqrtA^2;
% Compute the corrected mean motion (rad/s)
n = sqrt(MU/A^3) + delta_n;
% Compute the time that has passed since the time of ephemeris.
delta_time = time - toe;
% Compue the mean anomaly
M = M0 + n*delta_time;
% Compute the eccentric anomaly using Kepler's equation
E = M;
ratio = 1;
while abs(ratio) > 10^-8
E_error = E - e*sin(E) - M;
E_deriv = 1 - e*cos(E);
ratio = E_error/E_deriv;
E = E - ratio;
end
E = E + ratio;
% Compute the true anomaly.
nu = atan2(sqrt(1 - e^2)*sin(E) / (1 - e*cos(E)),...
(cos(E) - e) / (1 - e*cos(E)) );
% Compute the nominal argument of latitude.
phi = nu + omega;
% Compute the corrected argument of latitude.
arg_lat = phi + Cus*sin(2*phi) + Cuc*cos(2*phi);
% Compute the corrected radius.
radius = A*(1 - e*cos(E)) + Crs*sin(2*phi) ...
+ Crc*cos(2*phi);
% Compute the corrected inclination.
incline = i0 + i_dot*delta_time ...
+ Cis*sin(2*phi) + Cic*cos(2*phi);
% Compute the corrected longitude of the ascending node.
loan = OMEGA + (OMEGA_dot - EARTH_ROT_RATE)...
*delta_time - EARTH_ROT_RATE*toe;
% Compute the orbital in-plane x-y position (ECI)
x_eci = radius * cos(arg_lat);
y_eci = radius * sin(arg_lat);
% Compute the ECEF WGS84 coordinates
ecef(1, 1) = x_eci * cos(loan) - y_eci * cos(incline) * sin(loan);
ecef(1, 2) = x_eci * sin(loan) + y_eci * cos(incline) * cos(loan);
ecef(1, 3) = y_eci * sin(incline);
end
function ned = ecef2ned(ecef, ref_ecef, WGS)
% Converts the WGS84 Earth-Centered, Earth-Fixed (ECEF) Cartesian
% coordinates to the North-East-Down (NED) Cartesian coordinates using
% the reference position and latitude and longitude for the NED frame.
ref_lla = ecef2lla(ref_ecef, WGS);
delta_ecef = ecef - ref_ecef;
north = -sind(ref_lla(1)).*cosd(ref_lla(2)).*delta_ecef(1) ...
- sind(ref_lla(1)).*sind(ref_lla(2)).*delta_ecef(2) ...
+ cosd(ref_lla(1)).*delta_ecef(3);
east = -sind(ref_lla(2)).*delta_ecef(1) ...
+ cosd(ref_lla(2)).*delta_ecef(2);
down = -cosd(ref_lla(1)).*cosd(ref_lla(2)).*delta_ecef(1) ...
- cosd(ref_lla(1)).*sind(ref_lla(2)).*delta_ecef(2) ...
- sind(ref_lla(1)).*delta_ecef(3);
ned = [north, east, down];
end
  1 个评论
Michael Van de Graaff
Michael Van de Graaff 2022-3-1
That's waaaay to much to look though (points for not giving too little info though). how about you give like 3 rows of one of your nine 901x6 arrays
so i'm expecting something like
sample_data = [timestamp1,a1,b1,c1,d1,e1; timestamp2,a2,b2,c2,d2,e2; timestamp3,a3,b3,c3,d3,e3];
and somehting similar (without a,b,c,d,e) for the unique_time vector. (make sure there's at least one match) Then we can at least zero in on what's actually going on.
Also, you could extract the Second, Minute, Hour,Day, Month, Year, from the timestamps and just do a nested loop. Have you looked at the datetime dovumentation?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Gyan Vaibhav
Gyan Vaibhav 2023-11-18
Hi Terry,
I understand that you are trying to find the occurrences of your unique timestamps in the mentioned array. I assume that you are trying to extract the 2nd to 4th columns for the matched timestamps, into Data_mat.
The given code can be modified as follows to achieve the expected results:
%% Find each unique time step
time_mat = [mat1(:, 1), mat2(:, 1), mat3(:, 1), mat4(:, 1), mat5(:, 1),mat6(:, 1), mat7(:, 1), mat8(:, 1), mat9(:, 1)];
time_linear = time_mat(:);
time_sort = sort(time_linear);
time_unique = unique(time_sort);
% Initialize Data_mat as an empty matrix
Data_mat = [];
% Loop through each unique time stamp
for i = 1:length(time_unique)
% Initialize a temporary matrix to store rows matching the current time stamp
temp_mat = [];
% Loop through each data matrix
for j = 1:9
% Find rows with matching time stamp in the j-th matrix
matching_rows = mat{j}(mat{j}(:, 1) == time_unique(i), :);
% Append matching rows to the temporary matrix
temp_mat = [temp_mat; matching_rows(:, 2:end)];
end
% If there are matching rows, append them to Data_mat
% Otherwise, append [0, 0, 0] to Data_mat
if ~isempty(temp_mat)
Data_mat = [Data_mat; temp_mat];
else
Data_mat = [Data_mat; 0, 0, 0];
end
end
This approach gives us Data_mat that now contains the extracted rows of data based on the matching timestamps as expected.
Hope this helps.
Thanks
Gyan

类别

在 Help Center 和 File Exchange 中查找有关 Satellite Mission Analysis 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by