Analyze GPS Satellite Visibility

Open Live Script

This example shows how to simulate and analyze GPS satellite visibility at specified receiver positions and times. Use live script controls to set different parameters for the satellite simulation.

Specify Simulation Parameters

Specify the RINEX navigation message file, duration in hours, and time between samples in seconds of the simulation. Also specify the position of the receiver in geodetic coordinates and the mask angle, or minimum elevation angle, of the receiver.

rinexfile = "GODS00USA_R_20211750000_01D_GN.rnx";
numHours = 24;
dt =60; % s
latitude =42.3013162; % deg
longitude = -71.3782972; % deg
altitude = 50; % m
recPos = [latitude longitude altitude]; % [deg deg m]
maskAngle = 5; % deg

Use the rinexread function to obtain the simulation start time and the initial satellite orbital parameters from the RINEX file. Use only the first entry for each satellite in the RINEX file.

data = rinexread(rinexfile);
gpsData = data.GPS;
[~, idx] = unique(gpsData.SatelliteID);
navmsg = gpsData(idx,:);
startTime = navmsg.Time(1);

Generate Satellite Visibilities

Using the parameters, generate the satellite visibilities as a matrix of logical values. Each row in the matrix corresponds to a time step, and each column corresponds to a satellite. To plot visibilities, iterate through the time vector calculating the satellite positions and look angles based on GNSS constellation simulation.

secondsPerHour = 3600;
timeElapsed = 0:dt:(secondsPerHour*numHours);
t = startTime + seconds(timeElapsed);

satPos = gnssconstellation(t(1),navmsg);
numSats = size(satPos,1);
numSamples = numel(t);
az = zeros(numSamples,numSats);
el = zeros(numSamples,numSats);
vis = false(numSamples,numSats);
satIDs = 1:numSats;

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

for ii = 1:numel(t)
    satPos = gnssconstellation(t(ii),navmsg);
    [az(ii,:),el(ii,:),vis(ii,:)] = lookangles(recPos,satPos,maskAngle);
    set(sp,AzimuthData=az(ii,vis(ii,:)), ...
        ElevationData=el(ii,vis(ii,:)), ...
        LabelData=satIDs(vis(ii,:)));
    drawnow limitrate
end

Plot Results

Use the logical matrix to generate a satellite visibility chart and plot the total number of visible satellites at each time step. In general, at least four satellites must be visible to compute a positioning solution.

visPlotData = double(vis);
visPlotData(visPlotData == false) = NaN; % Hide invisible satellites.
visPlotData = visPlotData + (0:numSats-1); % Add space to satellites to be stacked.
colors = colororder;

figure
plot(t,visPlotData," .",Color=colors(1,:))
yticks(1:numSats)
yticklabels(string(satIDs))
grid on
ylabel("Satellite ID")
xlabel("Time")
title("Satellite Visibility Chart")
axis tight

numVis = sum(vis,2);
figure
area(t,numVis)
grid on
xlabel("Time")
ylabel("Number of satellites visible")
title("Number of GPS satellites visible")
axis tight