skyplot

Create a GNSS sensor model as a gnssSensor (Navigation Toolbox) System object™.

gnss = gnssSensor; 

Specify the position and velocity of the sensor. Simulate the sensor readings and get status from visible satellites. Store the azimuth and elevation angles as vectors.

pos = [0 0 0]; 
vel = [0 0 0]; 
[~, ~, status] = gnss(pos, vel); 
satAz = status.SatelliteAzimuth; 
satEl = status.SatelliteElevation; 

Create random local elevation masks, with a maximum elevation of 30 degrees, to act as the local environment.

rng(8)
terrainMaskElevations = 30*rand(1,12); % elevations (degrees)
terrainMaskEdges = [0 24 48 100 132 180 204 240 276 288 300 312 360]; % azimuth edges (degrees)

Plot the satellite positions with the elevation masks.

skyplot(satAz,satEl,MaskElevation=terrainMaskElevations,MaskAzimuthEdges=terrainMaskEdges); 

Animate the trajectory of satellite positions over time from a GNSS sensor.

Initialize the sky plot figure. Specify the relevant time-stepping information.

skyplotHandle = skyplot(0,0); 

numHours = 12; 
dt = 100; 
numSeconds = numHours * 60 * 60; 
numSimSteps = numSeconds/dt;

Create a GNSS sensor model as a gnssSensor (Navigation Toolbox) System Object™.

gnss = gnssSensor('SampleRate', 1/dt); 

Iterate through the time steps and do the following:

for i = 1:numSimSteps 

    [~, ~, status] = gnss([0 0 0],[0 0 0]); 

    satAz = status.SatelliteAzimuth; 
    satEl = status.SatelliteElevation; 

    set(skyplotHandle,'AzimuthData',satAz,'ElevationData',satEl);
    
    drawnow 
end 

Load the azimuth and elevation data from a logfile generated by an Adafruit® GPS satellite sensor. The data provided in this example contains the azimuth and elevation of each satellite and the pseudorandom noise (PRN) codes. Store these values as vectors.

load('gpsHWInfo','hwInfo') 
satAz = hwInfo.SatelliteAzimuths; 
satEl = hwInfo.SatelliteElevations; 
prn = hwInfo.SatellitePRNs; 

Separate the satellites based on the PRN codes. To correlate each position with a group, create a categorical array. For this set of satellites, only the ones with PRNs less than 32 are used in the positioning solution.

isUnused = (prn > 32);
group = categorical(isUnused,[false true],["Used in Positioning Solution" "Unused"]);

Visualize the satellites and specify the categorical groups in the GroupData name-value argument. Specify the PRN as the label for each point. Show the legend.

skyplot(satAz,satEl,prn,GroupData=group) 
legend('Used','Unused')

Specify the receiver position, RINEX navigation file, mask angle, time step size, and number of hours of data to sample from the RINEX file.

recPos = [42 -71 50];
navfile = "GODS00USA_R_20211750000_01D_GN.rnx";
maskAngle = 25;
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. Then, create the time vector t.

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

Initialize 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 satellite trajectories as an animation by creating a skyplot and updating the AzimuthData and ElevationData properties.

figure
sp = skyplot(allAz(1,:),allEl(1,:),satID,MaskElevation=maskAngle);
for idx = 1:size(allAz, 1)
    set(sp,AzimuthData=allAz(1:idx,:),ElevationData=allEl(1:idx,:));
    drawnow limitrate
end