Commit 0cd511e28c01f3da848b243d107fdd230ae3888a

Authored by Dany Dumont
1 parent 65cc596e
Exists in master and in 2 other branches mingan, simon

Segmentation des trajectoires

Showing 1 changed file with 141 additions and 28 deletions   Show diff stats
tr_map.m
... ... @@ -16,7 +16,6 @@ function tr_map
16 16 % device and the name the mat file. The number of devices must
17 17 % equal the number specified at the first line.
18 18 %
19   -%
20 19 % Example of infile:
21 20 %
22 21 % 3
... ... @@ -26,7 +25,6 @@ function tr_map
26 25 % spot17_20140912a
27 26 % spot19
28 27 % spot19_20140910a
29   -%
30 28 %
31 29 % Other m-files required: Image Processing Toolbox
32 30 % Subfunctions: none
... ... @@ -45,14 +43,14 @@ function tr_map
45 43 curdir = pwd;
46 44 matdir = [curdir,'/mat'];
47 45  
48   -if ~exist(matdir)
  46 +if ~exist(matdir,'dir')
49 47 disp(' ERROR : The /mat directory does not exist. Make')
50 48 disp(' sure mat files exist and that they ')
51 49 disp(' are placed in a directory named mat. ')
52 50 return
53 51 end
54 52  
55   -if exist('infile')
  53 +if exist('infile','file')
56 54 fid = fopen('infile');
57 55 else
58 56 disp(' ERROR : You need an input file named infile ')
... ... @@ -68,7 +66,7 @@ else
68 66 return
69 67 end
70 68  
71   -ndstr = fgetl(fid);
  69 +ndstr = fgetl(fid); % Read the 1st line of infile
72 70 nd = str2double(ndstr);
73 71  
74 72 lonmin = 360;
... ... @@ -81,21 +79,30 @@ for n = 1:nd
81 79 A(n).file = fgetl(fid);
82 80 load([matdir,'/',A(n).file]);
83 81  
  82 + clear sgmnts
  83 +
84 84 [~,ds] = legs(data(:,3),data(:,2)); % distance (nm)
85   - dt = data(2:end,1) - data(1:end-1,1); % time difference (days)
  85 + dt = diff(data(:,1)); % time difference (days)
  86 + dtmax = 65./(60.*24); % maximum time for valid speed
86 87 v = 1000.*nm2km(ds)./(dt.*86400); % average speed (m s-1)
  88 + vmask = dt < dtmax; % mask for speed (0,1)
  89 + %vmask(vmask == 0) = NaN;
87 90 travel = nm2km(sum(ds)); % km
  91 + sgmnts = find(dt > 1.1);
88 92  
89 93 A(n).data = data;
90 94 A(n).ds = ds;
91 95 A(n).dt = dt;
92 96 A(n).v = v;
  97 + A(n).vmask = vmask;
93 98 A(n).travel = travel;
  99 + A(n).sgmnts = sgmnts;
94 100  
95 101 lonmin = min(lonmin,min(A(n).data(:,2)));
96 102 lonmax = max(lonmax,max(A(n).data(:,2)));
97 103 latmin = min(latmin,min(A(n).data(:,3)));
98 104 latmax = max(latmax,max(A(n).data(:,3)));
  105 +
99 106 end
100 107  
101 108 %% Compute and display relevant quantities
... ... @@ -105,29 +112,50 @@ disp(&#39; &#39;)
105 112 for n = 1:nd
106 113 disp([' ',A(n).dev]);
107 114 disp(' ')
108   - disp([' Number of points : ',num2str(length(A(n).data(:,2)))]);
109   - disp([' Total travel : ',num2str(A(n).travel),' km']);
110   - disp([' Initial position : ',num2str(A(n).data(1,2)),'E, ',num2str(A(n).data(1,3)),'N']);
111   - disp([' Final position : ',num2str(A(n).data(end,2)),'E, ',num2str(A(n).data(end,3)),'N']);
112   - disp([' Deployment time : ',datestr(A(n).data(1,1))]);
113   - disp([' Beach time : ',datestr(A(n).data(end,1))]);
114   - disp([' Total time : ',num2str(A(n).data(end,1) - A(n).data(1,1)),' days']);
115   - disp([' Average speed : ',num2str(1000.*A(n).travel./(86400.*(A(n).data(end,1) - A(n).data(1,1)))),' m/s']);
  115 + disp([' Number of points : ',num2str(length(A(n).data(:,2)))]);
  116 + disp([' Total travel : ',num2str(A(n).travel),' km']);
  117 + disp([' Initial position : ',num2str(A(n).data(1,2)),'E, ',num2str(A(n).data(1,3)),'N']);
  118 + disp([' Final position : ',num2str(A(n).data(end,2)),'E, ',num2str(A(n).data(end,3)),'N']);
  119 + disp([' Deployment time : ',datestr(A(n).data(1,1)),' GMT']);
  120 + disp([' Beach time : ',datestr(A(n).data(end,1)),' GMT']);
  121 + disp([' Total time : ',num2str(A(n).data(end,1) - A(n).data(1,1)),' days']);
  122 + disp([' Average speed : ',num2str(1000.*A(n).travel./(86400.*(A(n).data(end,1) - A(n).data(1,1)))),' m/s']);
  123 + disp([' Maximum speed : ',num2str(max(A(n).v)),' m/s']);
  124 + disp([' Number of segments : ',num2str(length(A(n).sgmnts))]);
116 125 disp(' ')
117 126 end
118 127 disp(' -------------------------------------- ')
119 128  
120 129  
121   -%% Create the map
  130 +%% Make figures
122 131 % Track and land colors
123   -trk(1,:,1) = [0.8 0.2 0.2];
124   -trk(1,:,2) = [0.2 0.8 0.2];
125   -trk(1,:,3) = [0.2 0.2 0.8];
126   -trk(1,:,4) = [0.4 0.2 0.8];
127   -trk(1,:,5) = [0.2 0.4 0.8];
  132 +if nd <= 6
  133 + % Couleurs chaudes
  134 + trk(1,:,1) = [1.0 0.0 0.0];
  135 + trk(1,:,2) = [0.8 0.3 0.2];
  136 + trk(1,:,3) = [0.7 0.8 0.1];
  137 + % Couleurs froides
  138 + trk(1,:,4) = [0.4 0.2 0.8];
  139 + trk(1,:,5) = [0.2 0.4 0.8];
  140 + trk(1,:,6) = [0.0 0.0 0.8];
  141 +else
  142 + trk = nan.*ones(1,3,100);
  143 + trk(1,1,:) = 0.0; % R
  144 + trk(1,2,:) = 0.0; % G
  145 + trk(1,3,:) = 0.0; % B
  146 +end
128 147  
129 148 %land_color = [241 235 144]./255;
130   -land_color = [0.7 0.7 0.7];
  149 +land_color = [0.75 0.75 0.75];
  150 +
  151 +% GSL
  152 +% lonmin = -63;
  153 +% lonmax = -56;
  154 +% latmin = 46.8;
  155 +% latmax = 51.5;
  156 +
  157 +% Locations of interest
  158 +lon_oh = -60.394274; lat_oh = 48.051471; % Old Harry
131 159  
132 160 % Close-up
133 161 tol = (lonmax - lonmin)./10;
... ... @@ -136,23 +164,108 @@ lonmax = lonmax + tol;
136 164 latmin = latmin - tol;
137 165 latmax = latmax + tol;
138 166  
  167 +%% Map
139 168 figure(1)
140 169 m_proj('mercator', ...
141 170 'longitudes',[lonmin lonmax], ...
142 171 'latitudes',[latmin latmax]);
143 172  
144   -m_gshhs_f('patch',land_color);
  173 +%m_gshhs_f('patch',land_color);
  174 +m_gshhs_h('patch',land_color,'edgecolor',land_color);
145 175 hold on
  176 +%m_gshhs_h('color',land_color);
146 177  
147   -for n = 1:nd
148   - h(n) = m_plot(A(n).data(:,2),A(n).data(:,3),'-o','LineWidth',2, ...
149   - 'Color',trk(1,:,n), ...
150   - 'MarkerSize',1, ...
151   - 'MarkerFaceColor',trk(1,:,n));
  178 +for n = 1:nd
  179 + nsg = length(A(n).sgmnts);
  180 + if nsg > 0
  181 + ms = 1;
  182 + for m = 1:nsg
  183 + 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', ...
  184 + 'LineWidth',1, ...
  185 + 'Color',trk(1,:,n), ...
  186 + 'MarkerSize',3, ...
  187 + 'MarkerFaceColor',trk(1,:,n));
  188 + ms = A(n).sgmnts(m);
  189 + end
  190 + else
  191 + m_plot(A(n).data(:,2),A(n).data(:,3),'o', ...
  192 + 'LineWidth',1, ...
  193 + 'Color',trk(1,:,n), ...
  194 + 'MarkerSize',3, ...
  195 + 'MarkerFaceColor',trk(1,:,n));
  196 + end
152 197 end
153 198  
  199 +% Add Old Harry on the map
  200 +m_plot(lon_oh,lat_oh,'ok', ...
  201 + 'MarkerSize',8, ...
  202 + 'MarkerFaceColor',[0 0 0]);
  203 +
154 204 m_grid('box','fancy','linestyle','none','fontsize',13,'fontname','Verdana');
155 205 hold off
156 206  
157 207 print(gcf,'map.png','-dpng','-painters');
158   -print(gcf,'map.eps','-depsc2','-painters');
159 208 \ No newline at end of file
  209 +print(gcf,'map.eps','-depsc2','-painters');
  210 +
  211 +%% Drifter speed
  212 +figure(2)
  213 +subplot(2,1,1)
  214 +xmin = 1e16;
  215 +xmax = 0;
  216 +
  217 +for n = 1:nd
  218 + nsg = length(A(n).sgmnts);
  219 + if nsg > 0
  220 + ms = 1;
  221 + for m = 1:nsg
  222 + plot(A(n).data(ms+1:A(n).sgmnts(m)-1,1),A(n).v(ms+1:A(n).sgmnts(m)-1),'ok', ...
  223 + 'MarkerSize',8, ...
  224 + 'MarkerFaceColor',trk(1,:,n));
  225 + ms = A(n).sgmnts(m);
  226 + hold on
  227 + end
  228 + else
  229 + plot(A(n).data(2:end,1),A(n).v,'ok', ...
  230 + 'MarkerSize',8, ...
  231 + 'MarkerFaceColor',trk(1,:,n));
  232 + hold on
  233 +
  234 +
  235 + end
  236 + xmin = min(xmin,A(n).data( 1,1));
  237 + xmax = max(xmax,A(n).data(end,1));
  238 +end
  239 +hold off
  240 +datetick
  241 +set(gca,'FontSize',14,'PlotBoxAspectRatio',[1.0 0.3 1.0]);
  242 +xlim([xmin xmax])
  243 +ylabel('Buoy velocity (m/s)')
  244 +
  245 +%% Time interval
  246 +subplot(2,1,2)
  247 +for n = 1:nd
  248 + nsg = length(A(n).sgmnts);
  249 + if nsg > 0
  250 + ms = 1;
  251 + for m = 1:nsg
  252 + plot(A(n).data(ms+1:A(n).sgmnts(m)-1,1),A(n).dt(ms+1:A(n).sgmnts(m)-1).*1440,'ok', ...
  253 + 'MarkerSize',8, ...
  254 + 'MarkerFaceColor',trk(1,:,n));
  255 + ms = A(n).sgmnts(m);
  256 + hold on
  257 + end
  258 + else
  259 + plot(A(n).data(2:end,1),A(n).dt.*1440,'ok', ...
  260 + 'MarkerSize',8, ...
  261 + 'MarkerFaceColor',trk(1,:,n));
  262 + end
  263 +end
  264 +hold off
  265 +datetick
  266 +set(gca,'FontSize',14,'PlotBoxAspectRatio',[1.0 0.3 1.0]);
  267 +xlim([xmin xmax])
  268 +ylim([0 3000])
  269 +ylabel('Time interval (min)')
  270 +
  271 +print(gcf,'stats.png','-dpng','-painters');
  272 +print(gcf,'stats.eps','-depsc2','-painters');
160 273 \ No newline at end of file
... ...