How to get satellite track for multi constellations?

I have this code, running:

data = rinexread(file);

[~,satIdx] = unique(data.GPS.SatelliteID);

navmsg = data.GPS(satIdx,:);

startTime = navmsg.Time(1);

secondsPerHour = 3600;

dt = 60; % seconds

numHours = 4;

timeElapsed = 0:dt:(secondsPerHour*numHours);

t = startTime + seconds(timeElapsed);

maskAngle = 10;

%Generate satellite visibilities

allSatSys = ["GPS","Galileo","GLONASS","BeiDou","QZSS","SBAS"];

satLetter = ["G","E","R","C","J","S"];

sp = skyplot([],[],MaskElevation=maskAngle);

figure

hold on;

for ii = 1:numel(allSatSys)

satSys = allSatSys(ii);

% Get the satellite positions of the current satellite system.

[navmsg,satIDs] = exampleHelperParseGNSSFile(gnssFileType,satSys,file);

satPos = gnssconstellation(startTime,navmsg);

[az,el,vis] = lookangles(recPos,satPos,maskAngle);

% Combine the satellite system symbol letter with the satellite ID.

satIDLabel = arrayfun(@(x) sprintf("%c%02d",satLetter(ii),x),satIDs);

% Create a categorical array to associate the current values with the

% current satellite system in the skyplot.

satGroup = categorical(repmat(ii,numel(satIDLabel),1),1:numel(allSatSys),allSatSys);

% Update the skyplot.

set(sp, ...

AzimuthData=[sp.AzimuthData(:); az(vis)], ...

ElevationData=[sp.ElevationData(:); el(vis)], ...

LabelData=[sp.LabelData(:); satIDLabel(vis)], ...

GroupData=[sp.GroupData; satGroup(vis)]);

numSats = numel(navmsg.SatelliteID);

allAz = NaN(numel(t),numSats);

allEl = allAz;

for idx = 1:numel(t)

[satPos,~,satID] = gnssconstellation(t(idx),RINEXData=navmsg);

[az,el,vis] = lookangles(recPos,satPos,maskAngle);

allAz(idx,:) = az;

allEl(idx,:) = el;

end

legend

end

%Functions

function [navmsg,satIDs] = exampleHelperParseGNSSFile(gnssFileType,satSys,file)

switch gnssFileType

case "Rinex"

switch satSys

case "GPS"

navmsg = rinexread(file);

% For RINEX files, extract GPS data and use only the first entry for each satellite.

gpsData = navmsg.GPS;

[~,idx] = unique(gpsData.SatelliteID);

navmsg = gpsData(idx,:);

case "Galileo"

navmsg = rinexread(file);

% For RINEX files, extract Galileo data and use only the first entry for each satellite.

galData = navmsg.Galileo;

[~,idx] = unique(galData.SatelliteID);

navmsg = galData(idx,:);

case "GLONASS"

navmsg = rinexread(file);

% For RINEX files, extract GLONASS data and use only the first entry for each satellite.

gloData = navmsg.GLONASS;

[~,idx] = unique(gloData.SatelliteID);

navmsg = gloData(idx,:);

case "BeiDou"

navmsg = rinexread(file);

% For RINEX files, extract BeiDou data and use only the first entry for each satellite.

bdsData = navmsg.BeiDou;

[~,idx] = unique(bdsData.SatelliteID);

navmsg = bdsData(idx,:);

case "QZSS"

navmsg = rinexread(file);

% For RINEX files, extract QZSS data and use only the first entry for each satellite.

qzsData = navmsg.QZSS;

[~,idx] = unique(qzsData.SatelliteID);

navmsg = qzsData(idx,:);

case "SBAS"

file = "GOP600CZE_R_20211750000_01D_SN.rnx";

navmsg = rinexread(file);

% For RINEX files, extract SBAS data and use only the first entry for each satellite.

sbaData = navmsg.SBAS;

[~,idx] = unique(sbaData.SatelliteID);

navmsg = sbaData(idx,:);

end

satIDs = navmsg.SatelliteID;

end

end

And I get the skyplot for GPS only from the example. Below is the example of the code, it is working fine, but I cannot get two together.

maskAngle = 10;

dt = 60; % seconds

numHours = 4;

%Read the navigation file and get the GPS data of all satellites captured

%in the file

data = rinexread(navfile);

[~,satIdx] = unique(data.GPS.SatelliteID);

navmsg = data.GPS(satIdx,:);

%Set the starting time to the initial time of the navigation message, the

%create the T vector

startTime = navmsg.Time(1);

secondsPerHour = 3600;

timeElapsed = 0:dt:(secondsPerHour*numHours);

t = startTime + seconds(timeElapsed);

%Initialise vectors for azimuth and elevation. Then, collect azimuth and

%elevation data at times t for all satellites.

numSats = numel(navmsg.SatelliteID);

allAz = NaN(numel(t),numSats);

allEl = allAz;

for idx = 1:numel(t)

[satPos,~,satID] = gnssconstellation(t(idx),RINEXData=navmsg);

[az,el,vis] = lookangles(recPos,satPos,maskAngle);

allAz(idx,:) = az;

allEl(idx,:) = el;

end

%Mark all satellites below the horizon with an elevation less than 0 as

%missing

allEl(allEl < 0) = missing;

%Display the satellit trajectories as an animationby creating a skyplot and

%updating the AzimuthData and ElevationData properties

figure

h = skyplot(allAz(1,:),allEl(1,:),satID,MaskElevation=maskAngle);

for idx = 1:size(allAz, 1)

% h.Labels = ["GPS"]

set(h,AzimuthData=allAz(1:idx,:),ElevationData=allEl(1:idx,:));

drawnow limitrate

end

How to have a single plot of satellites track for all constellations?

[Kautuk Raj]

You can consider sharing the file being used in the statement data = rinexread(file) for better understanding and reproduction.

Answer 1

[Ronit]

Hello Joe,

To generate a single ‘skyplot’ that includes satellite tracks for all GNSS constellations, you need to ensure that you are processing and plotting data for all satellite systems together. You can use the following points to achieve this:

1. Combine Satellite Data: Process each GNSS constellation separately but store their azimuth and elevation data in a combined manner.
2. Update the Skyplot: Use a single ‘skyplot’ object and update it with the combined data for all satellite systems.

combinedAz = [];
combinedEl = [];
combinedLabels = [];

% Read the navigation file
data = rinexread(navfile);

% Create a skyplot figure
figure;
sp = skyplot([], [], MaskElevation=maskAngle);

% Process each satellite system
for ii = 1:numel(allSatSys)
    satSys = allSatSys(ii);
    [navmsg, satIDs] = exampleHelperParseGNSSFile("Rinex", satSys, navfile);
    
    startTime = navmsg.Time(1);
    t = startTime + seconds(timeElapsed);
    numSats = numel(navmsg.SatelliteID);
    allAz = NaN(numel(t), numSats);
    allEl = allAz;
    
    for idx = 1:numel(t)
        [satPos, ~, satID] = gnssconstellation(t(idx), RINEXData=navmsg);
        [az, el, vis] = lookangles(recPos, satPos, maskAngle);
        allAz(idx, :) = az;
        allEl(idx, :) = el;
    end
    
    % Append data for the current constellation
    combinedAz = [combinedAz; allAz(:)];
    combinedEl = [combinedEl; allEl(:)];
    satIDLabel = arrayfun(@(x) sprintf("%c%02d", satLetter(ii), x), satIDs, 'UniformOutput', false);
    combinedLabels = [combinedLabels; repmat(satIDLabel, numel(t), 1)];
end

% Update the skyplot with combined data
set(sp, ...
    'AzimuthData', combinedAz, ...
    'ElevationData', combinedEl, ...
    'LabelData', combinedLabels);

% Add legend
legend(sp, allSatSys);

Please refer to this MATLAB Answer for more details on plotting multiple sets using a single ‘skyplot’: https://www.mathworks.com/matlabcentral/answers/1610590-using-skyplot-and-hold-on?s_tid=answers_rc1-3_p3_MLT

I hope it helps you query!