GatherEddiesIntoAtlantic.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% What if moving all the oceanic mesoscale eddies into the Atlantic Ocean ?
% 
% This idea results from a discussion with one of my colleagues who did
% not believe oceanic mesoscale eddies cover nearly a third of the ocean 
% surface at any time. The data calculation such as summing up all the
% areas of eddies is effective but sometimes boring. Thus, how about 
% visualizing it by gathering them together in one place ? The vast Atlantic 
% Ocean covers about 29 percent of the Earth's water surface area according
% to Wikipedia (https://www.wikiwand.com/en/Atlantic_Ocean), making it the
% perfect domain for this kind of demonstration. 
%
% Data:
%   Global Mesoscale Eddy Track Atlas Product released by AVISO
%     https://www.aviso.altimetry.fr/en/data/products/value-added-products/global-mesoscale-eddy-trajectory-product.html
%
% Dependancy:
%   Lon/Lat coordinate data for boundary of the Atlantic Ocean, including islands.
%     Provided by http://www.marineregions.org/
%
%   The M_Map toolbox
%     and associated Lon/Lat coordinate data for coastlines provided by the
%       M_Map toolbox.
%     https://www.eoas.ubc.ca/~rich/map.html
%   land_or_ocean.m
%     https://ww2.mathworks.cn/matlabcentral/fileexchange/45268-land_or_ocean-m
%
% Published at:
%   https://github.com/chouj/MoveEddiesIntoAtlantic


% Load Global Mesoscale Eddy Track Atlas Product released by AVISO
% https://www.aviso.altimetry.fr/en/data/products/value-added-products/global-mesoscale-eddy-trajectory-product.html

% Juliandate
j1=ncread('tracks_AVISO_DT2014_daily_web.nc','j1');    

% Eddy polarity, 1 for cyclonic eddies and -1 for anticyclonic eddies
cyc=ncread('tracks_AVISO_DT2014_daily_web.nc','cyc');

% Eddy center longitudes 
lon=ncread('tracks_AVISO_DT2014_daily_web.nc','lon');

% Eddy center latitudes
lat=ncread('tracks_AVISO_DT2014_daily_web.nc','lat');

% Eddy radius (Km), eddy occupying area is represented as a circle having radius of L.
L=ncread('tracks_AVISO_DT2014_daily_web.nc','L');     

% for circle drawing
ang=[-pi:0.001:pi]; % Angles
x=cos(ang);y=sin(ang);

% Load coastline data from M_Map package.
% https://www.eoas.ubc.ca/~rich/map.html
% Modify the path according to yours.
load D:\MATLAB\toolbox\MatlabToolAdded\m_map\private\m_coasts.mat

% Load Lon/Lat coordinate data for the boundary of the Atlantic Ocean,
% including islands.Provided by http://www.marineregions.org/ and a .mat
% version generated by me can be downloaded here: 
% https://github.com/chouj/MoveEddiesIntoAtlantic/raw/master/AtlanticOceanBoundary.mat
load AtlanticOceanBoundary.mat

% Rough boundary of the Atlantic Ocean
bA=[-67.3,-60;20,-60;20,10;-6.75,10;-6.75,43.64;-1.56,43.64;-6.75,58.49;...
    -0.95,61.13;-37.51,70.97;-52.27,48.47;-76.71,42.85;-83,30.4;-112.8,30.4...
    ;-83,8.93;-75.5,8.93;-60,8;-60,-30;-67.3,-30;-67.3,-60];

% Convert Julian Date to Calendar days
dayj=min(j1):max(j1);
day=datetime(dayj,'convertfrom','juliandate');

% n indicates the indices for eddies detected on the first day of the
%   temporal range of dataset.
% Modify it according to your interest.
n=find(j1==min(j1));

% finding eddies whose center fall within the Atlantic Ocean
nn=inpolygon(lon(n)-360,lat(n),bA(:,1),bA(:,2));

% Figure plotting: draw all the eddy circles in the Atlantic Ocean
% Map Projection
m_proj('robinson','lon',[-100 290]);

for i=1:length(n)
    if cyc(n(i))==1&nn(i)
        m_plot((x*L(n(i))/110-360+lon(n(i))),(y*L(n(i))/110+lat(n(i))),'r'); % red for anticyclones
    elseif cyc(n(i))==-1&nn(i)
        m_plot((x*L(n(i))/110-360+lon(n(i))),(y*L(n(i))/110+lat(n(i))),'b'); % blue for cyclones
    end
    hold on
end

m_coast('patch',[0.7 0.7 0.7]);hold on
m_grid('linest','-');
% m_plot(bA(:,1),bA(:,2),'g-','linewidth',2);

% Generate property matrix for eddies in the Atlantic Ocean.
% A eddy list for both the original ones that not need to move and those 
% after moved.
EddyIn(:,2)=lat(n(nn));     % Lat
EddyIn(:,3)=L(n(nn))/110;   % Eddy radius in degree
EddyIn(:,1)=lon(n(nn))-360; % Lon
EddyIn(:,4)=cyc(n(nn));     % Eddy polarity

% Generate property matrix for eddies outside the Atlantic Ocean.
% Have to move them. Moving waiting list.
EddyOut(:,2)=lat(n(~nn));
EddyOut(:,3)=L(n(~nn))/110;
EddyOut(:,1)=lon(n(~nn))-360;
EddyOut(:,4)=cyc(n(~nn));

% Original eddy property matrix before moving eddies
EddyOutO=EddyOut;
EddyInO=EddyIn;

% Eddy moving starts

j=0; % count how many eddies have been moved.
while length(EddyOut(:,1))>0 % Continue if there still are eddies waiting for moving.
    
    % Randomly generate coordinate for eddy center after moved.
    Eddy.lon=rand(1)*[max(bA(:,1))-min(bA(:,1))]+min(bA(:,1));
    Eddy.lat=rand(1)*[max(bA(:,2))-min(bA(:,2))]+min(bA(:,2));
    
    % Diagnose whether the generated position fall within the region of the
    % Atlantic Ocean and whether it is on land or not
    if inpolygon(Eddy.lon,Eddy.lat,bA(:,1),bA(:,2))&land_or_ocean(Eddy.lat,Eddy.lon,10)
        clear disttoB disttoE disttoC 
        
        % Calculate the distances from generated coordinate to 1) the 
        % boundary of the Atlantic ocean
        disttoB=((Eddy.lon-bAd(:,1)).^2+(Eddy.lat-bAd(:,2)).^2).^0.5;
        
        % Calculate the distances from generated coordinate to 2) the 
        % coastline
        disttoC=((Eddy.lon-ncst(:,1)).^2+(Eddy.lat-ncst(:,2)).^2).^0.5;
        
        % Calculate the distances from generated coordinate to 3) all 
        % existing eddy circles (not eddy centers).
        disttoE=((Eddy.lon-EddyIn(:,1)).^2+(Eddy.lat-EddyIn(:,2)).^2).^0.5-EddyIn(:,3);
        
        % Continue if the generated coordinate does not fall within any 
        % existing eddy circles
        if min(disttoE)>0
            clear nnn sortL I
            
            % Find whether there are eddies in the moving waiting list that
            % can be fit into the minimum distance
            nnn=find(EddyOut(:,3)<=min([min(disttoB),min(disttoE),min(disttoC)]));
            if length(nnn)>0
                
                % Sort it by the radius. Find the largest eddy.
                [sortL,I]=sort(EddyOut(nnn,3),'descend');
                
                % Figure plotting
                if EddyOut(nnn(I(1)),4)==1 % Anticyclone
                    
                    % meganta circles for eddy before moving
                    m_plot(x*EddyOut(nnn(I(1)),3)+EddyOut(nnn(I(1)),1),y*EddyOut(nnn(I(1)),3)+EddyOut(nnn(I(1)),2),'m');
                    
                    % transparent red circles for eddy after moved
                    he=m_plot(x*EddyOut(nnn(I(1)),3)+Eddy.lon,y*EddyOut(nnn(I(1)),3)+Eddy.lat,'r');
                    he.Color(4) = 0.2;
                    
                elseif EddyOut(nnn(I(1)),4)==-1 % Cyclone
                    
                    % cyan circles for eddy before moving
                    m_plot(x*EddyOut(nnn(I(1)),3)+EddyOut(nnn(I(1)),1),y*EddyOut(nnn(I(1)),3)+EddyOut(nnn(I(1)),2),'c');
                    
                    % transparent blue circles for eddy after moved
                    he=m_plot(x*EddyOut(nnn(I(1)),3)+Eddy.lon,y*EddyOut(nnn(I(1)),3)+Eddy.lat,'b');
                    he.Color(4) = 0.2;
                end
                
                j=j+1; % Counting
             
                lEddyIn=length(EddyIn(:,1)); % how many eddies are in the Atlantic Ocean before moving this eddy.
                
                % Fill the information of the newly moved eddy into the list.
                EddyIn(lEddyIn+1,:)=EddyOut(nnn(I(1)),:);
                EddyIn(lEddyIn+1,1)=Eddy.lon;
                EddyIn(lEddyIn+1,2)=Eddy.lat;
                
                % Record its original information in a new matrix.
                % List them following the moving sequence so you know which
                % one is moved first and which one is the last.
                EddyOutS(j,:)=EddyOut(nnn(I(1)),:);
                
                % delete the information of the newly moved eddy in moving
                % waiting list.
                EddyOut(nnn(I(1)),:)=[];
                disp(length(EddyOut(:,1))) % display how many eddies have not been moved yet.
            end
        end
    end
end

% Figure plotting
title({'What If moving all oceanic mesoscale eddies into the vast domain of the Atlantic Ocean ?'...
        ['¡ª¡ªGlobal Mesoscale Eddy Trajectory Atlas Product (AVISO) ',datestr(day(1),'yyyy-mm-dd')]});

t1=m_text(90,-76.5,{'Anticyclones({\color{red}O}) and cyclones({\color{blue}O}) in the Atlantic Ocean.'...
    'Anticyclones({\color{magenta}O}) and cyclones({\color{cyan}O}) outside the Atlantic Ocean that have been moved'...
    'to new positions indicated by transparent red and blue circles in the Atlantic Sector.'});
t1.FontSize=9;
t1.BackgroundColor='w';
t1.HorizontalAlignment='center';

% export_fig EddiesMovedIntoAtlanticOcean.png -png -r600 -q111 -transparent


%%%%%%%%%%%%%%%%%%%%% Animation Part %%%%%%%%%%%%%%%%%%%%

filename = 'MovingAllEddiesIntoAtlantic1.gif';

% Prepare the new video file.
vidObj = VideoWriter('MovingAllEddiesIntoAtlantic1.mp4','MPEG-4');
vidObj.FrameRate=24;
vidObj.Quality=100;
open(vidObj);

m_proj('robinson','lon',[-100 290]);

% Eddy property matrix for those eddies after moved, for example, their new
% positions.
EddyInS=EddyIn(length(EddyInO)+1:end,:);

% 2-D grid points for 2-D interpolation. 
[X,Y]=meshgrid([1,101],1:length(EddyInS));
[Xi,Yi]=meshgrid(1:101,1:length(EddyInS));

% 2-D interpolation. To generate moving trajectories: positions during
% moving.
EddyMlon=interp2(X,Y,[EddyOutS(:,1),EddyInS(:,1)],Xi,Yi,'linear');
EddyMlat=interp2(X,Y,[EddyOutS(:,2),EddyInS(:,2)],Xi,Yi,'linear');

% back and forth 
EddyMlon=[EddyMlon,fliplr(EddyMlon(:,1:100))];
EddyMlat=[EddyMlat,fliplr(EddyMlat(:,1:100))];

% Animation producing
for j=1:201
    disp(j);
    h = figure('Visible', 'off');
    a = axes('Visible','off');
    h.PaperPositionMode = 'auto';
    h.PaperUnits = 'points';
    h.Units = 'points';
    h.Position = [6 43.2 1000.4 506]; % modify accordingly
    
    % Draw all eddies in the Atlantic Ocean that do not need to move.
    for i=1:length(n)
        if cyc(n(i))==1&nn(i)
            m_plot((x*L(n(i))/110-360+lon(n(i))),(y*L(n(i))/110+lat(n(i))),'r');
        elseif cyc(n(i))==-1&nn(i)
            m_plot((x*L(n(i))/110-360+lon(n(i))),(y*L(n(i))/110+lat(n(i))),'b');
        end
        hold on
    end
    
    % Draw all eddies being moved.
    for i=1:length(EddyInS)
        if EddyInS(i,4)==1
            m_plot(x*EddyOutS(i,3)+EddyMlon(i,j),y*EddyOutS(i,3)+EddyMlat(i,j),'m');
        elseif EddyInS(i,4)==-1
            m_plot(x*EddyOutS(i,3)+EddyMlon(i,j),y*EddyOutS(i,3)+EddyMlat(i,j),'c');
        end
        hold on
    end
    
    m_coast('patch',[0.7 0.7 0.7]);hold on
    m_grid('linest','-');
    
    title({'What If moving all oceanic mesoscale eddies into the vast domain of the Atlantic Ocean ?'...
        ['¡ª¡ªGlobal Mesoscale Eddy Trajectory Atlas Product (AVISO) ',datestr(day(1),'yyyy-mm-dd')]});
    t1=m_text(90,-76.5,{'Anticyclones({\color{red}O}) and cyclones({\color{blue}O}) in the Atlantic Ocean.'...
        'Anticyclones({\color{magenta}O}) and cyclones({\color{cyan}O}) outside the Atlantic Ocean being moved.'...
        'Source Code: https://github.com/chouj/MoveEddiesIntoAtlantic'});
    t1.FontSize=9;
    t1.BackgroundColor='w';
    t1.HorizontalAlignment='center';
    
    % Output
    frame = getframe(h);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    % Write to the GIF File and MP4 File.
    if j == 1
        imwrite(imind,cm,filename,'gif','DelayTime',1.5, 'Loopcount',inf);
        for ii=1:36
        writeVideo(vidObj,frame);
        end
    elseif j==101
        imwrite(imind,cm,filename,'gif','DelayTime',2.5,'WriteMode','append');
        for ii=1:60
        writeVideo(vidObj,frame);
        end
    else
        imwrite(imind,cm,filename,'gif','DelayTime',0.05,'WriteMode','append');
        writeVideo(vidObj,frame);
    end
    close
end
close(vidObj);