tr_map.m 7.76 KB
function tr_map
% TR_MAP - Display and save a map showing drifter trajectories specified
% in the input file.
%
% Syntax:  tr_map
% 
% Inputs:  The function reads an ascii input file called infile, which
%          must be written according to the following rules:
%          - The first line indicates the number of trajectories or GPS
%            devices to plot.
%          - The second line indicates the name of the device (e.g. spot15).
%          - The third line indicates the name of the mat file containing
%            the trajectory data. This file must be located in a directory
%            called mat/, which is located where infile is.
%          - Subsequent lines indicate, in the same order, the name of the
%            device and the name the mat file. The number of devices must 
%            equal the number specified at the first line.
%
% Example of infile:
%
%          3
%          spot15
%          spot15_20140910a
%          spot17
%          spot17_20140912a
%          spot19
%          spot19_20140910a
% 
% Other m-files required: Image Processing Toolbox
% Subfunctions: none
% MAT-files required: none
% 
% See also: TR_MERGE, GPX2MAT
% 
% Author: Dany Dumont
% UQAR/ISMER
% email: dany_dumont@uqar.ca
% Website: http://www.ismer.ca/dumont-dany
% September 2014
% ______________________________________________________________________


curdir = pwd;
matdir = [curdir,'/mat'];

if ~exist(matdir,'dir')
    disp(' ERROR : The /mat directory does not exist. Make')
    disp('         sure mat files exist and that they     ')
    disp('         are placed in a directory named mat.   ')
    return
end

if exist('infile','file')
    fid = fopen('infile');
else
    disp(' ERROR : You need an input file named infile    ')
    disp('         which must be written according to the ')
    disp('         following rules:                       ')
    disp('         3                % # of trajectories   ')
    disp('         spot15           % device 1            ')
    disp('         spot15_20140910a % mat file of device 1')
    disp('         spot17           % device 2            ')
    disp('         spot17_20140912a % mat file of device 2')
    disp('         spot19           % device 3            ')
    disp('         spot19_20140910a % mat file of device 3')
    return
end

ndstr = fgetl(fid);         % Read the 1st line of infile
nd    = str2double(ndstr);

lonmin =  360;
lonmax = -360;
latmin =  90;
latmax = -90;

for n = 1:nd
    A(n).dev  = fgetl(fid);
    A(n).file = fgetl(fid);
    load([matdir,'/',A(n).file]);
    
    clear sgmnts
    
    [~,ds] = legs(data(:,3),data(:,2));         % distance (nm)
    dt     = diff(data(:,1));                   % time difference (days)
    dtmax  = 65./(60.*24);                      % maximum time for valid speed
    v      = 1000.*nm2km(ds)./(dt.*86400);      % average speed (m s-1)
    vmask  = dt < dtmax;                        % mask for speed (0,1)
    %vmask(vmask == 0) = NaN;
    travel = nm2km(sum(ds));                    % km
    sgmnts = find(dt > 1.1);
    
    A(n).data   = data;
    A(n).ds     = ds;
    A(n).dt     = dt;
    A(n).v      = v;
    A(n).vmask  = vmask;
    A(n).travel = travel;
    A(n).sgmnts = sgmnts;
    
    lonmin = min(lonmin,min(A(n).data(:,2)));
    lonmax = max(lonmax,max(A(n).data(:,2)));
    latmin = min(latmin,min(A(n).data(:,3)));
    latmax = max(latmax,max(A(n).data(:,3)));
    
end

%% Compute and display relevant quantities

disp(' -------------- Summary -------------- ')
disp(' ')
for n = 1:nd
    disp([' ',A(n).dev]);
    disp(' ')
    disp(['   Number of points   : ',num2str(length(A(n).data(:,2)))]);
    disp(['   Total travel       : ',num2str(A(n).travel),' km']);
    disp(['   Initial position   : ',num2str(A(n).data(1,2)),'E, ',num2str(A(n).data(1,3)),'N']);
    disp(['   Final position     : ',num2str(A(n).data(end,2)),'E, ',num2str(A(n).data(end,3)),'N']);
    disp(['   Deployment time    : ',datestr(A(n).data(1,1)),' GMT']);
    disp(['   Beach time         : ',datestr(A(n).data(end,1)),' GMT']);
    disp(['   Total time         : ',num2str(A(n).data(end,1) - A(n).data(1,1)),' days']);
    disp(['   Average speed      : ',num2str(1000.*A(n).travel./(86400.*(A(n).data(end,1) - A(n).data(1,1)))),' m/s']);
    disp(['   Maximum speed      : ',num2str(max(A(n).v)),' m/s']);
    disp(['   Number of segments : ',num2str(length(A(n).sgmnts))]);
    disp(' ')
end
disp(' -------------------------------------- ')


%% Make figures
% Track and land colors
if nd <= 6
    % Couleurs chaudes
    trk(1,:,1) = [1.0 0.0 0.0];
    trk(1,:,2) = [0.8 0.3 0.2];
    trk(1,:,3) = [0.7 0.8 0.1];
    % Couleurs froides
    trk(1,:,4) = [0.4 0.2 0.8];
    trk(1,:,5) = [0.2 0.4 0.8];
    trk(1,:,6) = [0.0 0.0 0.8];
else
    trk = nan.*ones(1,3,100);
    trk(1,1,:) = 0.0; % R
    trk(1,2,:) = 0.0; % G
    trk(1,3,:) = 0.0; % B
end

%land_color = [241 235 144]./255;
land_color = [0.75 0.75 0.75];

% GSL
% lonmin = -63;
% lonmax = -56;
% latmin =  46.8;
% latmax =  51.5;

% Locations of interest
lon_oh = -60.394274; lat_oh =  48.051471;   % Old Harry

% Close-up
tol = (lonmax - lonmin)./10;
lonmin = lonmin - tol;
lonmax = lonmax + tol;
latmin = latmin - tol;
latmax = latmax + tol;

%% Map
figure(1)
m_proj('mercator', ...
    'longitudes',[lonmin lonmax], ...
    'latitudes',[latmin latmax]);

%m_gshhs_f('patch',land_color);
m_gshhs_h('patch',land_color,'edgecolor',land_color);
hold on
%m_gshhs_h('color',land_color);

for n = 1:nd 
    nsg = length(A(n).sgmnts);
    if nsg > 0
        ms = 1;
        for m = 1:nsg
            m_plot(A(n).data(ms+1:A(n).sgmnts(m)-1,2),A(n).data(ms+1:A(n).sgmnts(m)-1,3),'o', ...
                'LineWidth',1, ...
                'Color',trk(1,:,n), ...
                'MarkerSize',3, ...
                'MarkerFaceColor',trk(1,:,n));
            ms = A(n).sgmnts(m);
        end
    else
        m_plot(A(n).data(:,2),A(n).data(:,3),'o', ...
            'LineWidth',1, ...
            'Color',trk(1,:,n), ...
            'MarkerSize',3, ...
            'MarkerFaceColor',trk(1,:,n));
    end
end

% Add Old Harry on the map
m_plot(lon_oh,lat_oh,'ok', ...
    'MarkerSize',8, ...
    'MarkerFaceColor',[0 0 0]);

m_grid('box','fancy','linestyle','none','fontsize',13,'fontname','Verdana');
hold off

print(gcf,'map.png','-dpng','-painters');
print(gcf,'map.eps','-depsc2','-painters');

%% Drifter speed
figure(2)
subplot(2,1,1)
xmin = 1e16;
xmax = 0;

for n = 1:nd
    nsg = length(A(n).sgmnts);
    if nsg > 0
        ms = 1;
        for m = 1:nsg
            plot(A(n).data(ms+1:A(n).sgmnts(m)-1,1),A(n).v(ms+1:A(n).sgmnts(m)-1),'ok', ...
                'MarkerSize',8, ...
                'MarkerFaceColor',trk(1,:,n));
            ms = A(n).sgmnts(m);
            hold on
        end
    else
        plot(A(n).data(2:end,1),A(n).v,'ok', ...
            'MarkerSize',8, ...
            'MarkerFaceColor',trk(1,:,n));
        hold on
        
        
    end
    xmin = min(xmin,A(n).data(  1,1));
    xmax = max(xmax,A(n).data(end,1));
end
hold off
datetick
set(gca,'FontSize',14,'PlotBoxAspectRatio',[1.0 0.3 1.0]);
xlim([xmin xmax])
ylabel('Buoy velocity (m/s)')

%% Time interval
subplot(2,1,2)
for n = 1:nd
    nsg = length(A(n).sgmnts);
    if nsg > 0
        ms = 1;
        for m = 1:nsg
            plot(A(n).data(ms+1:A(n).sgmnts(m)-1,1),A(n).dt(ms+1:A(n).sgmnts(m)-1).*1440,'ok', ...
                'MarkerSize',8, ...
                'MarkerFaceColor',trk(1,:,n));
            ms = A(n).sgmnts(m);
            hold on
        end
    else
        plot(A(n).data(2:end,1),A(n).dt.*1440,'ok', ...
            'MarkerSize',8, ...
            'MarkerFaceColor',trk(1,:,n));
    end
end
hold off
datetick
set(gca,'FontSize',14,'PlotBoxAspectRatio',[1.0 0.3 1.0]);
xlim([xmin xmax])
ylim([0 3000])
ylabel('Time interval (min)')

print(gcf,'stats.png','-dpng','-painters');
print(gcf,'stats.eps','-depsc2','-painters');