Newer
Older
function [F] = erosanimation(variable,varargin)
% Visualize output of the EROS landscape evolution model as animations (LEM)
%
%
% The following function library is required, which can be downloaded
% from e.g. the MATLAB file exchange:
%
% TopoToolbox - A MATLAB program for the analysis of digital elevation
% models. (https://github.com/wschwanghart/topotoolbox)
%
%
% SYNTAX
%
% B = erosanimation(variable)
% B = erosanimation(variable,pn,pv,...)
%
%
% DESCRIPTION
%
% erosanimation creates movie frames of landscape evolution either in
% profile view, map view or in 3d.
%
%
% INPUT (required)
%
% variable variable of interest (string)
% 'sediment' Sediment thickness
% 'water' Water depth
% 'discharge' Water discharge
% 'qs' Unit-sediment flux
% 'downward' Flow orientation
% 'stress' Shear stress
% 'slope' Stream slope
% 'capacity' Stream capacity
% 'stock' Sediment stock
% 'hum' Water discharge on the topography
% 'rain' Sources (>0) and sinks (-1) of water and sediment
%
% 'profile' custom profile, second argument needs to be a DEM (GRIDobj)
% 'sprofile' stream long profile, second argument needs to be a DEM (GRIDobj)
%
% INPUT (optional)
%
% Parameter name/value pairs (pn,pv,...)
%
% 'mode' visualization mode (string) (default: 'movie2')
% 'movie2' 2d movie of variable
% 'movie3' 3d movie of topographic evolution
%
% 'viewdir' view geometry specified as 2-element vector of azimuth
% and elevation (default: [45,45])
% only apllies to mode 'movie3'
%
%
%
% OUTPUT
%
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
%
% EXAMPLE
%
% Run the example that comes with the Eros download and:
%
% 1. make an 2d-animation of sediment thickness and use the returned
% frames to construct an animated .gif
%
% eros_template.m
% B = erosanimation('sediment');
% frames2gif(B,'sediment.gif',0.1)
%
% 2. plot the average sediment thickness versus time
%
% B = erosanimation('sediment','mode','average');
%
% 3. make an 3d-animation of topography and return frames
%
% B = erosanmimation('topo','mode','movie3');
%
% REFERENCES:
%
% Davy, P., & Lague, D. (2009). Fluvial erosion/transport equation of land-
% scape evolution models revisited. Journal of Geophysical Research, 114,
% 116. https://doi.org/10.1029/2008JF001146.
%
% Davy, P., Croissant, T., & Lague, D. (2017). A precipiton method to cal-
% culate river hydrodynamics, with applications to flood prediction, land-
% scape evolution models, and braiding instabilities. Journal of
% Geophysical Research: Earth Surface, 122, 14911512.
% https://doi.org/10.1002/2016JF004156
%
%
% Author: Juergen Mey (juemey[at]uni-potsdam.de)
% Date: 28. May, 2020
p = inputParser;
expectedInput_variable = {'topo','water','sediment','flux','qs',...
'discharge','downward','stress','hum','slope','capacity','stock','sprofile','profile'};
addRequired(p,'variable',@(x) any(validatestring(x,expectedInput_variable)));
if strcmp(variable,'sprofile') || strcmp(variable,'profile')
addRequired(p,'dem',@(x)isa(x,'GRIDobj'));
default_flowmin = 100;
addParameter(p,'flowmin',default_flowmin,@isnumeric);
parse(p,variable,varargin{:});
dem = p.Results.dem;
flowmin = p.Results.flowmin;
else
default_mode = 'movie2';
expectedInput_mode = {'movie2','movie3'};
addParameter(p,'mode',default_mode,@(x) any(validatestring(x,expectedInput_mode)));
default_viewdir = [45,45];
addParameter(p,'viewdir',default_viewdir,@isnumeric);
parse(p,variable,varargin{:});
mode = p.Results.mode;
viewdir = p.Results.viewdir;
end
switch variable
case 'topo'
filetype = 'alt';
iylabel = 'Elevation (m)';
colors = 'landcolor';
case 'sediment'
filetype = 'sed';
iylabel = 'Sediment thickness (m)';
colors = 'jet';
case 'water'
filetype = 'water';
iylabel = 'Water depth (m)';
colors = 'flowcolor';
case 'capacity'
filetype = 'capacity';
iylabel = 'Capacity';
colors = 'jet';
case 'discharge'
filetype = 'discharge';
iylabel = 'Water discharge (m^3/s)';
colors = 'flowcolor';
case 'flux'
filetype = 'flux';
iylabel = 'Water discharge (m^3/s)';
colors = 'flowcolor';
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
case 'downward'
filetype = 'downward';
iylabel = 'Mean settling velocity (m/s)';
colors = 'parula';
case 'hum'
filetype = 'hum';
iylabel = 'Water discharge on topography (m^3/s)';
colors = 'flowcolor';
case 'qs'
filetype = 'qs';
iylabel = 'Sediment flux (m^3/s)';
colors = 'jet';
case 'slope'
filetype = 'slope';
iylabel = 'Slope';
colors = 'parula';
case 'stock'
filetype = 'stock';
iylabel = 'Sediment stock (m^3)';
colors = 'jet';
case 'stress'
filetype = 'stress';
iylabel = 'Shear stress (Pa)';
colors = 'jet';
end
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
if strcmp(variable,'sprofile')
FD = FLOWobj(dem,'preprocess','c');
S = STREAMobj(FD,flowacc(FD)>flowmin);
S2 = modify(S,'interactive','reachselect');
marker = round(linspace(1,length(S2.distance),10));
H = dir('*.alt');
W = dir('*.water');
S = dir('*.sed');
D = dir('*.flux');
[~,index] = sortrows({H.date}.');
H = H(index);
W = W(index);
S = S(index);
D = D(index);
for i = 1:length(D)
[z,~] = fopengrd(D(i).name);
z(z==0)=NaN;
B(:,:,i) = z;
end
sed = grd2GRIDobj(S(1).name,dem);
w = waitbar(1/length(H),['Collecting movie frames ... ']);
for i=1:length(H)
h=figure;
subplot(2,1,1)
water2 = grd2GRIDobj(W(i).name,dem);
dem2 = grd2GRIDobj(H(i).name,dem);
sed2 = grd2GRIDobj(S(i).name,dem);
plotdz(S2,dem,'color',[0.9 0.9 0.9]);hold on
plotdz(S2,dem-sed,'color',[0.7 0.7 0.7])
plotdz(S2,dem2+water2,'color',[0 0.61 1])
plotdz(S2,dem2,'color',[1 0.5 0.1]);hold on
plotdz(S2,dem2-sed2,'color','k')
xlim([0 S2.distance(end)])
yl(i,1:2) = ylim;
ylim([yl(1,1),yl(1,2)+20]);
title(['Time = ',num2str(i),''])
xlabel('Distance (m)')
ylabel('Elevation (m)');
legend({'Initial topo','Initial bedrock','Water','Sediment','Bedrock'},'Location','northwest')
legend('boxoff')
subplot(2,1,2)
flux = grd2GRIDobj(D(i).name,dem);
flux.Z(flux.Z==0)=NaN;
imageschs(dem2,flux,'colormap','flowcolor','caxis',[nanmin(B(:)),nanmax(B(:))],'colorbarylabel','Water discharge (m^3/s)');
hold on;
plot(S2,'k--')
% scatter(S2.x(marker),S2.y(marker),'k')
text(S2.x(marker),S2.y(marker),num2str(round(S2.distance(marker)/1000)),'FontWeight','bold')
x0=10;
y0=10;
width=2200;
height=1900/2;
set(gcf,'position',[x0,y0,width,height])
set(gcf,'Visible','off')
F(i) = getframe(gcf);
close all
waitbar(i/length(H))
end
close(w)
elseif strcmp(variable,'profile')
H = dir('*.alt');
W = dir('*.water');
S = dir('*.sed');
D = dir('*.flux');
[~,index] = sortrows({H.date}.');
H = H(index);
W = W(index);
S = S(index);
D = D(index);
for i = 1:length(D)
[z,~] = fopengrd(D(i).name);
z(z==0)=NaN;
B(:,:,i) = z;
end
sed = grd2GRIDobj(S(1).name,dem);
[d,z,x,y] = demprofile(dem);
[~,sz0] = demprofile(sed,[],x,y);
w = waitbar(1/length(H),['Collecting movie frames ... ']);
for i=1:length(H)
h=figure;
subplot(2,1,1)
water2 = grd2GRIDobj(W(i).name,dem);
dem2 = grd2GRIDobj(H(i).name,dem);
sed2 = grd2GRIDobj(S(i).name,dem);
[~,hz] = demprofile(dem2,[],x,y);
[~,wz] = demprofile(water2,[],x,y);
[~,sz] = demprofile(sed2,[],x,y);
plot(d,z,'color',[0.9 0.9 0.9]);hold on
plot(d,z-sz0,'color',[0.7 0.7 0.7])
plot(d,hz+wz,'color',[0 0.61 1])
plot(d,hz,'color',[1 0.5 0.1]);hold on
plot(d,hz-sz,'color','k')
xlim([0 d(end)])
yl(i,1:2) = ylim;
ylim([yl(1,1),yl(1,2)+20]);
title(['Time = ',num2str(i),''])
xlabel('Distance (m)')
ylabel('Elevation (m)');
subplot(2,1,2)
flux = grd2GRIDobj(D(i).name,dem);
flux.Z(flux.Z==0)=NaN;
imageschs(dem2,flux,'colormap','flowcolor','caxis',[nanmin(B(:)),nanmax(B(:))],'colorbarylabel','Water discharge (m^3/s)');
hold on;
plot(x,y,'k--')
text(x(1)+10,y(1)+10,'left');
text(x(end)+10,y(end)+10,'right');
x0=10;
y0=10;
width=2200;
height=1900/2;
set(gcf,'position',[x0,y0,width,height])
set(gcf,'Visible','off')
F(i) = getframe(gcf);
close all
waitbar(i/length(H))
end
close(w)
else
T = dir('*.ini');
Z = dir(['*.',filetype]);
[t,~] = fread_timeVec(T.name,length(Z));
if isempty(t)
t=1:length(Z);
end
if isnan(t)
t=1:length(Z);
end
[~,index] = sortrows({Z.date}.');
Z = Z(index);
for i = 1:length(Z)
[z,~] = fopengrd(Z(i).name);
z(z==0)=NaN;
B(:,:,i) = z;
end
switch mode
case 'movie2'
H = dir('*.alt');
Z = dir(['*.',filetype]);
[~,index] = sortrows({H.date}.');
H = H(index);
Z = Z(index);
w = waitbar(1/length(H),['Collecting movie frames ... ']);
for i = 1:length(H)-1
h = grd2GRIDobj(H(i+1).name);
z = grd2GRIDobj(Z(i+1).name);
z.Z(z.Z==0)=NaN;
%imageschs(h,z,'colormap',colors,'caxis',[0,1],'colorbarylabel',iylabel);
imageschs(h,z,'colormap',colors,'caxis',[nanmin(B(:)),nanmax(B(:))],'colorbarylabel',iylabel);
%set(gca,'ColorScale','log')
title(['Time = ',num2str(t(i)),''])
x0=10;
y0=10;
width=2200;
height=1900;
set(gcf,'position',[x0,y0,width,height])
set(gcf,'Visible','off')
F(i) = getframe(gcf);
close all
waitbar(i/length(H))
end
close(w)
case 'movie3'
H = dir('*.alt');
[~,index] = sortrows({H.date}.');
H = H(index);
w = waitbar(1/length(H),['Collecting movie frames ... ']);
for i = 1:length(H)
h = grd2GRIDobj(H(i).name);
[xm,ym] = getcoordinates(h);
axis off
surface(xm,ym,h.Z,'EdgeColor','none');colorbar
view(viewdir(1),viewdir(2))
axis equal
c = colorbar;
c.Label.String = 'Elevation (m)';
colormap(landcolor)
caxis([nanmin(B(:)),nanmax(B(:))])
title(['Time = ',num2str(t(i)),''])
x0=10;
y0=10;
width=2200;
height=1900;
set(gcf,'position',[x0,y0,width,height])
set(gcf,'Visible','off')
F(i) = getframe(gcf);
close all
waitbar(i/length(H))
end
close(w)
end