How can I plot the interpolation of a set of points whose x coordinates range from [-pi, pi] without having the line to go back and forth the figure (which has the axis from - pi to pi)?

Basically I'm trying to do the ground track of a satellite and I want to plot in real time (using drawnow) the interpolation of a set of calculated points representing the projection of the satellite's position in space on Earth. The longitude axis goes from -180 to 180 degrees. If I make matlab simply run:

for i = 2 : n-1                  % n is the number of points
plot(longitude(i-1: i+1), latitude(i-1: i+1), 'r'); hold on; drawnow
end

the plot will be covered with almost horizontal lines since the satellite will reach a point close to longitude = 180 and then it will reappear close to longitude = -180, on the other side of the figure. With the code I wrote before matlab will connect the former point with the latter with a line that goes all accross the figure from right to left. I've actually solved this problem in most cases writing a series of "if"s but I'd like to know if there is a way to avoid this problem without complicating the code with all of those "if"s. I'm asking this because of course i'd like to avoid simply plotting the points like this:

for i = 2 : n-1                  % n is the number of points
plot(longitude(i), latitude(i), 'r.'); hold on; drawnow
end

Yes, it woud solve the problem but it's uglier :D

[Image Analyst]

Can you make it easy for us to help you by uploading them in a .mat file?

save('Answers.mat', 'longitude', 'latitude');

or else attach your original data file and code to read it in.

Then use the paper clip icon to upload them.

[Riccardo Basile]

function [] = GroundTrackGFX(latitude, glongitude, date, Nt, inc)
gtf = figure('position', [25 200 800 500], 'NumberTitle', 'off',...
    'Name', 'Ground track');
M = imread('Earth texture.jpg');
S = imread('Satellite icon.png');
imagesc('XData', [-180 180],'YData', [90 -90],'CData', M); 
daspect([1 1 1]); hold on
xlabel('Longitude [deg]'); ylabel('Latitude [deg]');
axis([-180 180 -90 90]); box on;
ax = gca; ax.TickDir = 'out';

% Assuming that the satellite is always launched by the longitude for which
% the X axis of the IRF passes, the first position of the satellite 
% depends on the launch date selected by the user.

plot(glongitude(2), latitude(2), 'g.', 'markersize', 16); hold on
title(date(2));
I = imagesc('XData', [glongitude(2)-3 glongitude(2)+3], 'YData',...
    [latitude(2) - 3 latitude(2) + 3], 'CData', S);
if (inc == 90)                       % This is the polar orbit case, done by simply plotting the points.
    for i = 3 : Nt-1
        title(date(i));
        delete(I);
        I = imagesc('XData', [glongitude(i)-3 glongitude(i)+3],...
            'YData', [latitude(i) - 3 latitude(i) + 3], 'CData', S);
        plot(glongitude(i), latitude(i), 'r.'); hold on
        drawnow
    end
else
     for i = 3 : Nt-1
        title(date(i));
        delete(I);
        I = imagesc('XData', [glongitude(i)-3 glongitude(i)+3],...
            'YData', [latitude(i) - 3 latitude(i) + 3], 'CData', S);
        if (glongitude(i) > 0 && glongitude(i+1) < 0)
            lat = [latitude(i-1) latitude(i) (latitude(i+1)+latitude(i))/2];
            lon = [glongitude(i-1) glongitude(i) 180];
            plot(lon, lat, 'r');
            hold on;
        end
        if (glongitude(i-1) > 0 && glongitude(i) < 0)
            lat = [(latitude(i-1)+latitude(i))/2 latitude(i) latitude(i+1)];
            lon = [-180 glongitude(i) glongitude(i+1)];
            plot(lon, lat, 'r');
            hold on;
        end
        if (glongitude(i-1) < 0 && glongitude(i+1) < 0)
            plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
            hold on;
        end
        if (glongitude(i-1) > 0 && glongitude(i+1) > 0)
            plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
            hold on;
        end
        if (glongitude(i) == 180)
            plot(glongitude(i-2: i), latitude(i-2: i), 'r');
            hold on;
        end
        if (glongitude(i) == -180)
            plot(glongitude(i: i+2), latitude(i: i+2), 'r');
            hold on;
        end
%         if (glongitude(i) < 0 && glongitude(i+1) > 0)
%             plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
%             hold on;
%         end
%         if (glongitude(i) == 0)
%             plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
%             hold on;
%         end
%         if (glongitude(i) < 0 && glongitude(i+1) > 0)
%             lat = [latitude(i-1) latitude(i) (latitude(i+1)+latitude(i))/2];
%             lon = [glongitude(i-1) glongitude(i) -180];
%             plot(lon, lat, 'r');
%             hold on
%         end
%         if (glongitude(i-1) < 0 && glongitude(i) > 0)
%             lat = [(latitude(i-1)+latitude(i))/2 latitude(i) latitude(i+1)];
%             lon = [180 glongitude(i) glongitude(i+1)];
%             plot(lon, lat, 'r');
%             hold on
%         end
%         if (glongitude(i) < 0 && glongitude(i+1) > 0)
%             plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
%             hold on
%         end
%         if (glongitude(i) > 0 && glongitude(i+1) < 0)
%             plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
%             hold on
%         end
        
%         if (latitude(i) < 0 && latitude(i+1) > 0)
%             plot(glongitude(i), latitude(i), 'r.');
%             hold on
%         end
%         if (latitude(i) > 0 && latitude(i+1) < 0)
%             plot(glongitude(i), latitude(i), 'r.');
%             hold on
%         end
%         if (latitude(i-1) < 0 && latitude(i) > 0)
%             plot(glongitude(i), latitude(i), 'r.');
%             hold on
%         end
%         if (latitude(i-1) > 0 && latitude(i) < 0)
%             plot(glongitude(i), latitude(i), 'r.');
%             hold on
%         end
        drawnow
    end
end
plot(glongitude(Nt), latitude(Nt), 'm.', 'markersize', 16); delete(I);
title(date(Nt));
dcm = datacursormode(gtf);
set(dcm,'DisplayStyle','window','Enable','on');
set(dcm,'UpdateFcn',  @(s,e) tipcallback(s,e, glongitude, date));
% Now the datacursormode update function is implemented:
    function output_txt = tipcallback(~,event_obj, glongitude, date)
    pos = get(event_obj,'Position');
    [lia, glonloc] = ismembertol(pos(1), glongitude, eps(pos(1)));
    if ((pos(1) == 180 && lia == 0) || (pos(1) == -180 && lia == 0))
        output_txt = "This point is only graphical. It's not part of the calculated points";
    else
    thedate = date(glonloc);                                
    
    output_txt{1} = ['Long: ' num2str(pos(1)),'   Lat: ' num2str(pos(2))];
    output_txt{2} = thedate;
    end
    end
end

The longitude vector I attached in answers.mat is actually the glongitude vector you see in the code (it stands for ground longitude and it's different from the longitude of the satellite in the inertial reference frame (but this is not important). Also the latitude vector has one more component than the longitude vector because I defined another longitude coordinate separately (but this is also not important, you can simply see longitude as x and latitude as y components of the points. You can see what I was trying to do with the "if"s I was talking about. Just to let you know... inc is the inclination of the orbit, Nt is the number of points and date is a string vector. The last function is also not important for what I asked. Thanks in advance!

[Riccardo Basile]

The problem is even better seen with this other set of points, of a more complex orbit. You can see the lines and the empty spots around longitude = 0. Also, I forgot to mention that the first point is not fisical and doen't have to be plotted!!

By the way I haven't uploaded all the code because this function is part of an app and it would be too messy and also usless to upload it.

I've attached the result I get with this set of points. If you remove those "if"s you'll get those lines every time the satellite passes from right to left and vice versa

Answer 1

[TADA]

how about something like that:

figure(1);
hold on;
currIdx = 1;
for i = [find(abs(diff(longitude')) > 10), length(longitude)] 
    idx = currIdx:i;
    if length(idx) < 2
        style = '.r';
    else
        style = 'r';
    end
    plot(longitude(idx), latitude(idx), style);
    currIdx = i+1;
end

It finds the indices where theres a big jump in the longitude, and cuts the plot there

[TADA]

the styling condition can be removed if you don't care about displaying the start of the satellite trajectory from the first point, as soon as there are at least 2 points in each segment, it'll plot normally

[Riccardo Basile]

First of all, thanks for the answer! But what does that "a." mean? Using the code as it is gives me error because it doesn't recognize it. However I eventually kind of solved the problem by covering every possible case with if structures (in the code below there's every case a part from the polar orbit case)

for i = 3 : Nt-1
        title(date(i));
        delete(I);
        I = imagesc('XData', [glongitude(i)-3 glongitude(i)+3],...
            'YData', [latitude(i) - 3 latitude(i) + 3], 'CData', S);
        if (glongitude(i) > 90 && glongitude(i+1) < 0)                         % from -180 to 180 (in)
            lat = [latitude(i-1) latitude(i) (latitude(i+1)+latitude(i))/2];
            lon = [glongitude(i-1) glongitude(i) 180];
            plot(lon, lat, 'r');
            hold on;
        end
        if (glongitude(i-1) > 0 && glongitude(i) < -90)                        % from -180 to 180 (out)
            lat = [(latitude(i-1)+latitude(i))/2 latitude(i) latitude(i+1)];
            lon = [-180 glongitude(i) glongitude(i+1)];
            plot(lon, lat, 'r');
            hold on;
        end
        if (glongitude(i-1) < 0 && glongitude(i) > 0 && glongitude(i) > 90)    % from 180 to -180 (out)
            lat = [(latitude(i-1)+latitude(i))/2 latitude(i) latitude(i+1)];
            lon = [180 glongitude(i) glongitude(i+1)];
            plot(lon, lat, 'r');
            hold on;
        end
        if (glongitude(i) < 0 && glongitude(i+1) > 0 && glongitude(i) < -90)   % from 180 to -180 (in)
            lat = [latitude(i-1) latitude(i) ((latitude(i+1)+latitude(i))/2)];
            lon = [glongitude(i-1) glongitude(i) -180];
            plot(lon, lat, 'r');
            hold on;
        end
        if (glongitude(i-1) < 0 && glongitude(i+1) < 0)                        % intern points
            plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
            hold on;
        end
        if (glongitude(i-1) > 0 && glongitude(i+1) > 0)                        % intern points
            plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
            hold on;
        end
        if (glongitude(i) == 180)                                              % case glongitude(i) = 180
            plot(glongitude(i-2: i), latitude(i-2: i), 'r');
            hold on;
        end
        if (glongitude(i) == -180)                                             % case glongitude(i) = -180
            plot(glongitude(i: i+2), latitude(i: i+2), 'r');
            hold on;
        end
        if (glongitude(i) == 0)                                                % case glongitude = 0
            plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
            hold on;
        end
        if (glongitude(i) < 0 && glongitude(i+1) > 0 && glongitude(i) > -90)   % from 0- to 0+
            plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
            hold on;
        end
        if (glongitude(i) > 0 && glongitude(i+1) < 0 && glongitude(i) < 90)    % from 0+ to 0-
            plot(glongitude(i-1: i+1), latitude(i-1: i+1), 'r');
            hold on;
        end
        drawnow
    end
end

[TADA]

sorry, I forgot to remove it

I loaded the data in the mat file you uploaded into a struct:

a = load('answers2.mat');

I now updated the answer and removed it

[TADA]

I see that in your code you disregard the single point, so you can use the code in this answer without the styling condition:

figure(1);
hold on;
currIdx = 1;
for i = [find(abs(diff(longitude')) > 10), length(longitude)] 
    idx = currIdx:i;
    plot(longitude(idx), latitude(idx), 'r');
    currIdx = i+1;
end

In your code you plot only a few data points at a time in a loop

this is much slower than plotting the whole data simultaniously because you have a lot of display updates. You wrote that this is a realtime system, so maybe it doesn't matter to you if you only plot one data poit at a time. However, if you refresh the whole figure each time, it's much faster to plot all of it together.

[Riccardo Basile]

It says the dimensions of concatenated arrays are not consistent. The diff(longitude) vector is shorter than longitude apparently (1 element shorter)

>> Untitled2
Error using horzcat
Dimensions of arrays being concatenated are not consistent.

Error in Untitled2 (line 7)
for i = [find(abs(diff(longitude)) > 10), a] % a = length(longitude)

Ah the problem seems to be the find command. I'm trying to write it in another way.

[TADA]

that's because longitude is a column vector. you should index using row vectors,

that's the reason I transposed the vector into a row vector:

longitude'

then when I concatenated horizontally the last index in longitude vector it worked:

indices = [find(abs(diff(longitude')) > 10), length(longitude)] 

indices =

           1         939        1630        1917

whereas if i try concating the indices vertically:

indices = [find(abs(diff(longitude)) > 10); length(longitude)] 

indices =

           1
         939
        1630
        1917
        

Thats all good, but the problem starts when you loop over the vector:

% i takes a row vector:
iteration = 1;
for i = [find(abs(diff(longitude')) > 10) length(longitude)] 
    istr = ['[' char(strjoin(string(num2str(i)), ';')) ']' ];
    disp(['iteration #' num2str(iteration) ' i = ' istr]);
    iteration = iteration + 1;
end

this is the output:

iteration #1 i = [1]
iteration #2 i = [939]
iteration #3 i = [1630]
iteration #4 i = [1917]

% i takes a column vector:
iteration = 1;
for i = [find(abs(diff(longitude)) > 10); length(longitude)] 
    istr = ['[' char(strjoin(string(num2str(i)), ';')) ']' ];
    disp(['iteration #' num2str(iteration) ' i = ' istr]);
    iteration = iteration + 1;
end

thats the output now:

iteration #1 i = [   1; 939;1630;1917]

bottom line, you need to use a row vector in the for loop, so you need to transpose that fist part somewhere:

for i = [find(abs(diff(longitude')) > 10) length(longitude)]
% or
for i = [find(abs(diff(longitude)') > 10) length(longitude)]
% or
for i = [find(abs(diff(longitude))' > 10) length(longitude)]
% or
for i = [find(abs(diff(longitude)) > 10)' length(longitude)]
    

either of those would work...

[Riccardo Basile]

Ok I see now. I don't think I can use this method to have an animation of the ground track but you gave me another idea... Thanks a lot!

[TADA]

Ok, I didn't realise you wanted an animation. so slow refresh is exactly what you wanted.

In that case, you can combine the two approaches:

figure(1);
ax = axes();
xlim(ax, [-180 180]);
ylim(ax, [-90 90]);
hold on;
discontinuousPoints = find(abs(diff(longitude)) > 10);
Nt = length(longitude);
for i = 1 : Nt-1
    if any(discontinuousPoints == i)
        continue;
    end
    plot(longitude([i i+1]), latitude([i i+1]), 'r');
    drawnow();
end

[Riccardo Basile]

It works great and it's more compact than my solution but still, there are blank spots near 180 and -180 degrees (I need the plot to be continous). I think I can't get around some of those "if"s.