I want to color code a satellite trajectory plot using Matlab

Question

I am creating a plot in Matlab that maps the azimuth and elevation of several GPS satellites over time. Below is an example of an image of the plotted GPS satellite trajectories.

However, I also want to plot the S_4 index (scintillation effect caused by the ionosphere) on the same plot. I want to color code each of the GPS satellite paths to show a high S_4 index as red, and a low S_4 index as blue. I have data for all of this, but I cannot figure out a way to color code each GPS path.

Currently I know the position of each satellite at every second of a day (this is how I am able to plot the trajectory). I also have S_4 index values at different times and for different satellites (the S_4 index matrix is not the same size as the satellite position matrix).

Any ideas?

This is the function I used. It uses as input a matrix called plotMat:

plotMat=[svID(satellite id number), gpsSec, elevation(radians), azimuth(radians)]

My function:

% find out from 'plotMat' if plotting satellite locations or trajectories in
% addition determine how many satellites are being tracked and how many
% samples for each satellite (# samples / satellite must always be equal)
gpsTime = plotMat(1,2);
i = 1;
t = gpsTime;
while ((i ~= size(plotMat,1)) & (t == gpsTime))
  i = i + 1;
  t = plotMat(i,2);
end
if (t == gpsTime)
  sats = i;
else
  sats = i - 1;
end;
samples = size(plotMat,1) / sats;
SVs = plotMat(1:sats,1);
elevation = plotMat(:,3);
azimuth = plotMat(:,4);
% initialize polar - plotting area
figure(10);clf;
axis([-1.4 1.4 -1.1 1.1]);
axis('off');
axis(axis);
hold on;
% plot circular axis and labels
th = 0:pi/50:2*pi;
x = [ cos(th) .67.*cos(th) .33.*cos(th) ];
y = [ sin(th) .67.*sin(th) .33.*sin(th) ];
plot(x,y,'color','w');
text(1.1,0,'90','horizontalalignment','center');
text(0,1.1,'0','horizontalalignment','center');
text(-1.1,0,'270','horizontalalignment','center');
text(0,-1.1,'180','horizontalalignment','center'); 
% plot spoke axis and labels
th = (1:6)*2*pi/12;
x = [ -cos(th); cos(th) ];
y = [ -sin(th); sin(th) ];
plot(x,y,'color','w');
text(-.46,.93,'0','horizontalalignment','center');
text(-.30,.66,'30','horizontalalignment','center');
text(-.13,.36,'60','horizontalalignment','center');
text(.04,.07,'90','horizontalalignment','center');
% plot titles
if (samples == 1)
  title('Satellite Position Plot');
  subtitle = sprintf('GPS time : %.2f sec',plotMat(1,2));
else
  title('Satellite Trajectory Plot');
  subtitle = sprintf('GPS time range : %.2f sec to %.2f sec', ...
                     plotMat(1,2),plotMat(size(plotMat,1),2));
  text(-1.6,1,'SVID/Last Position','color','r');
  text(-1.6,.9,'Positive Elevation','color','y');
  text(-1.6,.8,'Negative Elevation','color','b');
end
text(0,-1.3,subtitle,'horizontalalignment','center');
% plot trajectories (or positions) and label the last postion with the
% satellite SV id; in addition, last postions are in red, otherwise position
% elevations are yellow and negative elavations are blue
%
% loop through samples
for s = 1:samples       
  % plot each satellite location for that sample
  for sv = 1:sats
    % check if positive or negative elevation
    if (elevation((s - 1) * sats + sv) < 0)
      elNeg = 1;
    else
      elNeg = 0;
    end
    % convert to plottable cartesian coordinates
    el = elevation((s - 1) * sats + sv);
    az = azimuth((s - 1) * sats + sv);
    x = (pi/2-abs(el))/(pi/2).*cos(az-pi/2);
    y = -1*(pi/2-abs(el))/(pi/2).*sin(az-pi/2);
    % check for final sample
    if (s == samples)
      plot(x,y,'r*');
      text(x,y+.07,int2str(SVs(sv)), ...
           'horizontalalignment', ...
           'center','color','r');
    else
      % check for +/- elevation
      if (elNeg == 0)
        plot(x,y,'y.');
      else
        plot(x,y,'b.');
      end
    end
  end
end

axis equal;

On the previous code I haven't added the S_4 index. The following code compares the GPS seconds on the scintillation (S_4) data with the GPS seconds with the satellite position data. The matrix scint.mat refers to the scintillation data, while txinfo.mat refers to the satellite position data. This code creates a combination code that agrees on the GPS second and satellite identification number.

load('scint_324_first')
scint = scint';

load('txinfo_324_first_wrw')
txinfo = txinfo';

k = 0;

for i = 1:length(scint)
    disp(num2str(i))
    for j = 1:length(txinfo)
        if scint(i,2) == txinfo(j,2) && scint(i,15) == txinfo(j,8)
            k = k+1;
            combo(k,:) = [scint(i,:) txinfo(j,:)];
        end
    end
end

Answer 1

If the S_4 values change over time, shall the colour coding change within single lines?

It's maybe not the most efficient way, but you could draw just lines for each piece of the orbits. Like

s4code = [1 0 0; 0 0 1]; % blue & red 

for i = time_vector
 if s4(i) > something
  s = 2
 else 
  s = 1
 end
 line(x(i:(i+1)), y(i:(i+1)), 'color', s4code(s, :))
end

About orbit data and s4 data not having the same size ... make them have the same size ... interp1, for example ;-)

Answer 2

If you are plotting using a polar plot, try instead 3D Polar Plot (FEX code by Ken Garrard)

If you are using plot, try Colored line or scatter plot or cline.m (again FEX files)

the last plot was done by:

x = 1:300
color_line(x,sin(x/20),sin(x/20))