Analyze GPS Satellite Visibility

Open Live Script

This example shows how to simulate and analyze GPS satellite visibility at specified receiver positions and times using a GPS ephemeris or almanac file. Use the live script controls to set various parameters for the satellite simulation.

Specify Simulation Parameters

Specify the navigation message file type, start time, 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.

gnssFileType = "RINEX";
startTime = datetime(2021,06,24,04,00,00);
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

Get Satellite Orbital Parameters

Use the exampleHelperParseGNSSFile helper function to obtain the initial satellite orbital parameters and satellite IDs from a GNSS file. Because ephemeris data is valid for two hours before and after the time of ephemeris, you can use it to determine the exact positions of the satellites in orbit. In contrast, while almanac data, is valid for a period of 90 days, you can only use it for estimating the approximate positions of the satellites.

[navmsg,satIDs] = exampleHelperParseGNSSFile(gnssFileType);

The files used in this example are tied to the specified start time startTime. To use different RINEX files, use the RINEX files included in the Read Navigation and Observation Data from RINEX File example. For examples of how to download and use different SEM and YUMA almanac files, see the Read GPS Navigation Message Data from SEM Almanac File and Read GPS Navigation Message Data from YUMA Almanac File, respectively.

Once you have new ephemeris data or an almanac file, update the file variable of the exampleHelperParseGNSSFile helper function for the case that corresponds to your file type.

Generate Satellite Visibilities

Using your simulation 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 while calculating the satellite positions and look angles based on the GNSS constellation simulation.

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

numSats = numel(satIDs);
numSamples = numel(t);
az = zeros(numSamples,numSats);
el = zeros(numSamples,numSats);
vis = false(numSamples,numSats);

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

for ii = 1:numel(t)
    satPos = gnssconstellation(t(ii),navmsg,GNSSFileType=gnssFileType);
    [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

Helper Function

Use the exampleHelperParseGNSSFile helper function to obtain the initial satellite orbital parameters and satellite IDs from the GNSS file, based on the file type.

function [navmsg,satIDs] = exampleHelperParseGNSSFile(gnssFileType)
    switch gnssFileType
        case "RINEX"
            file = "GODS00USA_R_20211750000_01D_GN.rnx";
            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,:);
            satIDs = navmsg.SatelliteID;
        case "SEM"
            file = "semalmanac_2021-6-22.al3";
            navmsg = semread(file);
            satIDs = navmsg.PRNNumber;
        case "YUMA"
            file = "yumaalmanac_2021-6-22.alm";
            navmsg = yumaread(file);
            satIDs = navmsg.PRN;
    end
end