erosanimation.m

function [B,varargout] = erosanimation(variable,varargin)
% Visualize output of the EROS landscape evolution model (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 shows timeseries data either as grid average or as 2d/3d movie
%
%
% INPUT (required)
%
%   variable       variable of interest (string)
%                  'topo'  Topographic elevation
%                  '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
%
% INPUT (optional)
%
%   Parameter name/value pairs (pn,pv,...)
%
%   'mode'         visualization mode (string) (default: 'movie2')
%                  'average'    shows the evolution of the spatial average
%                               of the variable defined as required input
%                  '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
%
%   B              movie frames captured with modes 'movie2' and 'movie3'
%   B              (mode='average') 3d array of the variable of interest.
%
%
% OUTPUT (optional)
%
%   meanB       spatial average of variable through time           
%
% 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','qs','flux',...
    'discharge','downward','stress','hum','slope','capacity','stock'};
addRequired(p,'variable',@(x) any(validatestring(x,expectedInput_variable)));

default_mode = 'movie2';
expectedInput_mode = {'average','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;

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';
    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

% determine timesteps

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);
    B(:,:,i) = z;
    meanB(i)=mean(z(:));
end
switch mode
    case 'average'
                plot(t,meanB)
                ylabel(iylabel);
                xlabel('time');
                grid on
                B=meanB;
    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',[min(B(:)),200],'colorbarylabel',iylabel);
            imageschs(h,z,'colormap',colors,'caxis',[min(B(:)),max(B(:))],'colorbarylabel',iylabel);
            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)
        B = F;
    case 'movie3'
        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)
            h = grd2GRIDobj(H(i).name);
            %             z = grd2GRIDobj(Z(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)),''])
            set(gcf,'Visible','off')
            F(i) = getframe(gcf);
            close all
            waitbar(i/length(H))
        end
        close(w)
        f = figure;
        movie(f,F,1,10)
        close(f)
        B = F;
end