few updates

@FLOWobj/flipdir.m

@@ -0,0 +1,56 @@
+function FD = flipdir(FD)
+
+%FLIPDIR Flip direction of flow
+%
+% Syntax
+%
+% FDm = flipdir(FD)
+%
+% Description
+%
+% flipdir changes the direction of flow in a single or multi FLOWobj so
+% that flow moves upstream. In upstream direction, single FLOWobjs
+% become divergent and thus flipdir always returns a FLOWobj with
+% multiple flow directions.
+%
+% Input arguments
+%
+% FD FLOWobj
+%
+% Output arguments
+%
+% FDm flipped FLOWobj (multi)
+%
+% Example
+%
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% FD = FLOWobj(fillsinks(DEM),'multi');
+% [FD,S] = multi2single(FD,'minarea',500);
+% FDf = flipdir(FD);
+% H = ~STREAMobj2GRIDobj(S);
+% A = flowacc(FDf,H);
+% imageschs(DEM,min(A,1000))
+%
+%
+% See also: FLOWobj, FLOWobj/multi2single
+%
+% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
+% Date: 5. September, 2019
+
+
+FD.type = 'multi';
+ix = FD.ix;
+FD.ix = FD.ixc;
+FD.ixc = ix;
+
+% flip ordering
+FD.ix = FD.ix(end:-1:1);
+FD.ixc = FD.ixc(end:-1:1);
+if isempty(FD.fraction)
+ FD.fraction = ones(size(FD.ix));
+else
+FD.fraction = FD.fraction(end:-1:1);
+end
+FD = multi_normalize(FD);
+
+

@FLOWobj/flowvec.m

@@ -0,0 +1,52 @@
+function [U,V] = flowvec(FD,L)
+
+%FLOWVEC velocity vectors from FLOWobj
+%
+% Syntax
+%
+% [U,V] = flowvec(FD)
+% [U,V] = flowvec(FD,L)
+% 
+% Description
+%
+% flowvec returns the velocity components of the flow network in FD.
+% The components are normalized so that vector length is one. The
+% vector
+
+
+U = GRIDobj(FD);
+[X,Y] = getcoordinates(U);
+X = single(X./FD.cellsize);
+Y = single(Y./FD.cellsize);
+[X,Y] = meshgrid(X,Y);
+
+DX = X(FD.ixc) - X(FD.ix);
+DY = Y(FD.ixc) - Y(FD.ix);
+
+clear X Y
+
+if ~ismulti(FD)
+ %% Single flow directions
+ V = U;
+ U.Z(FD.ix) = DX;
+ V.Z(FD.ix) = DY;
+
+else
+ %% Multiple flow directions 
+ U.Z = reshape(accumarray(FD.ix,DX,[prod(U.size) 1],...
+ @mean,nan('single')),U.size);
+ V = U;
+ V.Z = reshape(accumarray(FD.ix,DY,[prod(U.size) 1],...
+ @mean,nan('single')),U.size); 
+end
+
+% Normalize vectors to have unit length
+LL = sqrt(U.^2 + V.^2);
+U = U./LL;
+V = V./LL;
+
+if nargin == 2
+ U = U.*L;
+ V = V.*L;
+end
+

@FLOWobj/multi2single.m

@@ -40,14 +40,16 @@
 %
 % Example
 % 
-% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
-% FD = FLOWobj(DEM,'preprocess','carve','mex',true);
-% DEM = imposemin(FD,DEM,0.0001);
-% S = STREAMobj(FD,'minarea',1000);
-% DEMc = DEM;
-% DEMc.Z(S.IXgrid) = DEMc.Z(S.IXgrid)-100; 
-% FD = FLOWobj(DEMc,'multi');
-% FD = multi2single(FD,'minarea',100);
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% FD = FLOWobj(DEM,'preprocess','carve','mex',true);
+% DEM = imposemin(FD,DEM,0.0001);
+% FD = FLOWobj(DEM,'multi');
+% [FD,S] = multi2single(FD,'minarea',100);
+% A = flowacc(FD);
+% imageschs(DEM,log(A),'colormap',flowcolor)
+% hold on
+% plot(S,'k')
+% hold off
 %
 % See also: FLOWobj
 %
@@ -60,6 +62,8 @@
 addParameter(p,'unit','pixels',@(x) ischar(validatestring(x, ...
  {'pixels', 'mapunits'})));
 addParameter(p,'channelheads',[])
+addParameter(p,'probability',false)
+addParameter(p,'randomize',false)
 
 parse(p,varargin{:});
 
@@ -70,6 +74,10 @@
 
  otherwise
 
+ if p.Results.randomize
+ FD = randomize(FD);
+ end
+
  if p.Results.minarea > 0 || ~isempty(p.Results.channelheads)
 
  if isempty(p.Results.channelheads)
@@ -117,8 +125,40 @@
 
  RR = (1:numel(FD.ix))';
  IX = double(FD.ix);
- S = sparse(RR,IX,FD.fraction,max(RR),max(IX));
- [~,ii] = max(S,[],1);
+ if ~p.Results.probability
+ S = sparse(RR,IX,FD.fraction,max(RR),max(IX));
+ [~,ii] = max(S,[],1);
+ else
+ [IXX,ix] = sort(IX);
+ c = FD.fraction(ix);
+
+ % cumulative probabilities for each edge leaving each giver
+ for r = 2:numel(IXX)
+ if IXX(r) == IXX(r-1)
+ c(r) = c(r)+c(r-1);
+ end
+ end
+
+ [~,a,b] = unique(IXX);
+ R = rand(size(a));
+ R = R(b);
+
+ II = c > R;
+ I = II;
+ for r = 2:numel(IXX)
+ if IXX(r) == IXX(r-1) && II(r-1)
+ I(r) = false;
+ end
+ end
+
+ frac = false(size(FD.fraction));
+ frac(ix) = I;
+ S = sparse(RR,IX,frac,max(RR),max(IX));
+ [~,ii] = max(S,[],1);
+
+
+
+ end
  I = false(size(RR));
  I(ii) = true;
 
@@ -161,4 +201,7 @@
 
 
 end
-end
+end
+
+
+

@FLOWobj/randomize.m

@@ -0,0 +1,63 @@
+function FD = randomize(FD)
+
+%RANDOMIZE randomize multiple flow directions
+%
+% Syntax
+%
+% FDr = randomize(FDm)
+%
+% Description
+%
+% randomize takes an instance of a multiple flowdirections FLOWobj and
+% randomizes the flow directions. This removes the slope dependency of
+% the fraction that each cell transfers to its downstream neighbors. 
+%
+% Input arguments
+%
+% FDm FLOWobj (multiple flow directions)
+%
+% Output arguments
+%
+% FDr FLOWobj (randomized multiple flow directions
+%
+% Example
+%
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% DEM = filter(DEM,'mean',[21 21]);
+% DEM = fillsinks(DEM);
+% 
+% ax(1) = subplot(1,2,1);
+% FD = FLOWobj(DEM,'multi');
+% A = flowacc(FD);
+% imageschs(DEM,log(A),'colormap',flowcolor)
+% 
+% ax(2) = subplot(1,2,2);
+% FDr = randomize(FD);
+% Ar = flowacc(FDr);
+% imageschs(DEM,log(Ar),'colormap',flowcolor)
+% 
+% linkaxes(ax,'xy')
+% % Now zoom in to explore differences in flow patterns
+%
+%
+% See also: FLOWobj, FLOWobj/multi2single, FLOWobj/flowconvergence
+%
+% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
+% Date: 23. May, 2019
+
+% If FD has single flow directions, do nothing
+if ~ismulti(FD)
+ return
+end
+
+% randomize the flow fraction using a cont. uniform distribution
+FD.fraction = rand(size(FD.fraction));
+
+% normalize so that the sum in each fraction is one
+FD = multi_normalize(FD);
+
+% Normalize code
+% [~,~,ix] = unique(FD.ix);
+% s = accumarray(ix,FD.fraction);
+% s = s(ix);
+% FD.fraction = FD.fraction./s;

@GRIDobj/curvature.m

@@ -29,6 +29,8 @@
 % (see function blockproc)
 % 'useparallel' true or {false}: use parallel computing toolbox
 % 'blocksize' blocksize for blockproc (default: 5000)
+% 'meanfilt' true or {false}: if true, preprocess DEM with 
+% [3x3] mean filter. 
 %
 % Output arguments
 %
@@ -74,8 +76,12 @@
 addParamValue(p,'useblockproc',false,@(x) isscalar(x));
 addParamValue(p,'blocksize',5000,@(x) isscalar(x));
 addParamValue(p,'useparallel',false,@(x) isscalar(x));
+addParamValue(p,'meanfilt',false,@(x) isscalar(x));
 parse(p,varargin{:});
 
+if p.Results.meanfilt
+ DEM = filter(DEM);
+end
 
 % create a copy of the DEM instance
 C = DEM;

@STREAMobj/cumsum.m

@@ -0,0 +1,73 @@
+function c = cumsum(S,c,direction)
+
+%CUMSUM cumulative sum on stream network
+%
+% Syntax
+%
+% c = cumsum(S,a)
+% c = cumsum(S,a,direction)
+%
+% Description
+%
+% cumsum calculates the cumulative sum along a stream network. By
+% default, cumsum accumulates the values in downstream direction.
+%
+% Input arguments
+%
+% S STREAMobj
+% a node-attribute list
+% direction 'downstream' (default) or 'upstream'
+%
+% Output
+%
+% c node attribute list (double)
+%
+% Example
+%
+% % Two step calculation of upslope area
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% % flow direction
+% FD = FLOWobj(DEM,'preprocess','carve','mex',true);
+% % stream network
+% S = STREAMobj(FD,'minarea',500);
+% a1 = hillslopearea(S,FD);
+% a2 = cumsum(S,a1);
+%
+% See also: STREAMobj/cumtrapz
+%
+% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
+% Date: 19. August 2019
+
+narginchk(2,3)
+
+if isa(c,'GRIDobj')
+ c = getnal(S,c);
+else
+ if ~isnal(S,c)
+ error('The second input argument must be a node attribute list.')
+ end
+end
+
+c = double(c);
+
+if nargin <= 2
+ direction = 'downstream';
+else 
+ direction = validatestring(direction,{'upstream','downstream'},'STREAMobj/cumsum','direction',3);
+end
+
+% get sorted edge list 
+ix = S.ix;
+ixc = S.ixc;
+
+% traverse stream network 
+switch direction
+ case 'upstream'
+ for r = numel(S.ix):-1:1
+ c(ix(r)) = c(ixc(r)) + c(ix(r));
+ end
+ otherwise
+ for r = 1:numel(S.ix)
+ c(ixc(r)) = c(ix(r)) + c(ixc(r));
+ end
+end

@STREAMobj/getnal.m

@@ -7,6 +7,7 @@
 % z = getnal(S,DEM)
 % [z1,z2,...] = getnal(S,DEM1,DEM2,...)
 % nalstruct = getnal(S,DEM1,DEM2,...,'struct')
+% z = getnal(S)
 %
 % Description
 %
@@ -21,6 +22,9 @@
 % to nan-punctuated vectors that can be plotted by custom functions,
 % use the function STREAMobj2XY.
 %
+% With only S as input argument, getnal returns a node-attribute list
+% of zeros (double). 
+%
 % Input
 %
 % S STREAMobj