some updates

- added GRIDobj/clip
- renamed FLOWobj2gradient to gradient (overloaded FLOWobj method)
- removed bug in STREAMobj/orientation

@FLOWobj/FLOWobj.m

@@ -567,7 +567,7 @@
  end
 
 
-  end
+ end
 
  else
  %% Multiple flow direction

@FLOWobj/dbentropy.m

@@ -4,18 +4,25 @@
 %
 % Syntax
 %
-% [H,S] = dbentropy(M,siz,ix)
+% H = dbentropy(FD);
+% [H,S] = dbentropy(FD);
+% [H,S] = dbentropy(FD,IX)
 %
 % Description
 %
-% dbentropy calculates the Shannon Entropy (H) of a digital elevation
-% model. H quantifies the uncertainty associated with each pixel in the
-% DEM to drain towards a specific outlet.
+% dbentropy calculates the drainage basin entropy H (Schwanghart and
+% Heckmann 2012) based on multiple flow directions stored in FLOWobj
+% FD. Drainage basin entropy is derived for each pixel and quantifies
+% the uncertainty to assign an outlet pixel. Pixel that drain towards a
+% single outlet pixel have low entropy values, whereas pixels with high
+% values could be assigned to two or more outlet pixels. Drainage basin
+% entropy relies on a probabilistic interpretation of multiple flow 
+% directions.
 %
 % Input parameters
 %
 % FD multiple flow direction FLOWobj
-% IX linear indices into the DEM from which FD was derived.
+% IX linear indices of outlets
 %
 % Output arguments
 % 
@@ -44,40 +51,42 @@
 % See also: drainagebasins, FLOWobj
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
-% Date: 07. October, 2016
+% Date: 21. March, 2019
 
 M = FLOWobj2M(FD);
 siz = FD.size;
 
-if nargin == 1;
+if nargin == 1
  ix = find(((sum(M,2)) == 0) & (sum(M,1)' > 0));
  nrix = numel(ix);
  % check the number of outlets and break the problems into chunks if
  % necessary
 
  chunksize = 100;
 
- if nrix > chunksize;
+ if nrix > chunksize
  chunks = true;
  else
  chunks = false;
  end
 
  % In addition, suggest to output only some of the larger fuzzy drainage basins 
- if nargout >= 2 && chunks;
-
-
+ if nargout >= 2 && chunks
+
  mess = ['The number of probability matrices is quite large (k = ' num2str(nrix) ').\n' ...
  'It is suggested to return only the probability matrices of\n'...
  'the outlets with the largest drainage basins. Please provide\n'...
  'a number of how many you want. The default is n = ' num2str(chunksize) '. n : '];
+
  pout = input(mess);
+
  if isempty(pout)
  pout = 10;
  end
 
- A = (speye(prod(siz))-M')\ones(prod(siz),1);
- [~,ixsort] = sort(A(ix),'descend');
+ % calculate flow accumulation
+ A = flowacc(FD);
+ [~,ixsort] = sort(A.Z(ix),'descend');
  ix = ix(ixsort);
 
  end
@@ -103,7 +112,7 @@
 % B is the membership (partition) matrix (?) (fuzzy partition matrix). 
 % B contains the membership grades. (calculated by membership function below)
 
-if nargout >= 2;
+if nargout >= 2
  Pc = struct('P',{},'IX',{},'CA50',{},'CA99',{},...
  'CA95',{},'CA66',{},'CA33',{},...
  'CA05',{},'CA01',{});
@@ -120,8 +129,9 @@
 
  if nargout >= 2
 
- for r = 1:nrix;
- Pc(r).P = full(reshape(P(:,r),siz));
+ for r = 1:nrix
+ Pc(r).P = GRIDobj(FD);
+ Pc(r).P.Z = full(reshape(P(:,r),siz));
  Pc(r).IX = ix(r);
  ptemp = sort(full(P(:,r)),'descend');
  [~,m] = unique(ptemp);
@@ -152,8 +162,12 @@
  H = zeros(siz);
  counter = 1;
  sP = zeros(prod(siz),1);
+
+ h = waitbar(0,'Starting');
+ iters = ceil(nrix/chunksize);
 
- for r = 1:ceil(nrix/chunksize);
+ for r = 1:iters
+ waitbar(r/iters,h,[num2str(r) ' / ' num2str(iters)]);
  cols = ((r-1)*chunksize + 1):r*chunksize;
  cols(cols>nrix) = [];
 
@@ -163,10 +177,11 @@
 
  if nargout >= 2
 
- for rr = cols(1):cols(end);
+ for rr = cols(1):cols(end)
 
  if rr<=pout
- Pc(rr).P = full(reshape(P(:,counter),siz));
+ Pc(rr).P = GRIDobj(FD);
+ Pc(rr).P.Z = full(reshape(P(:,counter),siz));
  ptemp = sort(Pc(rr).P(:),'descend');
  else
  Pc(rr).P = [];
@@ -187,7 +202,7 @@
 
 
  % reset counter if larger than chunksize
- if counter >= chunksize;
+ if counter >= chunksize
  counter = 1;
  else
  counter = counter+1;
@@ -197,6 +212,7 @@
  end
  end
  end 
+ close(h);
 % H = H + reshape(1 - sP.*log2(sP),siz);
 end
 

@FLOWobj/FLOWobj2gradient.m → @FLOWobj/gradient.m

@@ -1,10 +1,10 @@
-function G = FLOWobj2gradient(FD,DEM)
+function G = gradient(FD,DEM)
 
-%FLOWOBJ2GRADIENT gradient along flow direction
+%GRADIENT gradient along flow direction
 %
 % Syntax 
 %
-% G = FLOWobj2gradient(FD,DEM)
+% G = gradient(FD,DEM)
 %
 % Description
 %
@@ -29,12 +29,13 @@
 %
 % DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
 % FD = FLOWobj(DEM,'preprocess','c');
-% G = FLOWobj2gradient(FD,DEM);
+% G = gradient(FD,DEM);
 % imageschs(DEM,G)
 %
+% See also: GRIDobj/gradient8, GRIDobj/arcslope, STREAMobj/gradient
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
-% Date: 23. February, 2013
+% Date: 27. March, 2019
 
 validatealignment(FD,DEM)
 d = getdistance(FD.ix,FD.ixc,FD.size,FD.cellsize);

@FLOWobj/private/copy2GRIDobj.m

@@ -8,6 +8,6 @@
 pg = properties(G);
 pf = properties(FD);
 p = intersect(pg,pf);
-for r = 1:numel(p);
+for r = 1:numel(p)
  G.(p{r}) = FD.(p{r});
 end

@GRIDobj/clip.m

@@ -0,0 +1,68 @@
+function [DEM,C] = clip(DEM,C)
+
+%CLIP clip a GRIDobj with a polygon or another GRIDobj
+%
+% Syntax
+%
+% DEMc = clip(DEM,MS)
+% DEMc = clip(DEM,C)
+% [DEMc,Cc] = clip(...)
+% 
+% Description
+%
+% This function clips a GRIDobj based on a polygon feature MS or
+% another GRIDobj C. C should contain a logical array. If the
+% underlying class is numeric, DEM is clipped to the non-nan pixels in
+% C. If the underlying class is integer, then DEM will be clipped to
+% the nonzero pixels in C. Pixels in DEM will be nan if they are not
+% within the clip features. If DEM contains a logical or integer
+% matrix, then pixels outside the clip features will be zero.
+%
+% Input arguments
+%
+% DEM GRIDobj
+% MS mapping structure array (see polygon2GRIDobj)
+% C GRIDobj. If C does not spatially align with DEM, then C will
+% be resampled (see GRIDobj/resample) using nearest neighbor
+% interpolation
+% 
+% Output arguments
+%
+% DEMc clipped GRIDobj
+% C GRIDobj that can be used as clip feature 
+%
+% See also: GRIDobj/crop, GRIDobj/polygon2GRIDobj
+% 
+% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
+% Date: 27. Marchs, 2019
+
+
+% Process input
+if isstruct(C)
+
+ C = polygon2GRIDobj(DEM,C);
+elseif isa(C,'GRIDobj')
+
+ % Check the data format and apply different ways to get a logical array
+ if isfloat(C.Z)
+ C.Z = isnan(C.Z);
+ elseif isinteger(C.Z)
+ C.Z = C.Z == 0;
+ end
+
+ tf = validatealignment(DEM,C);
+ if ~tf
+ C = resample(C,DEM,'nearest');
+ end
+end
+
+% Process output
+if isfloat(DEM.Z)
+ DEM.Z(~C.Z) = nan;
+elseif isinteger(DEM.Z)
+ DEM.Z(~C.Z) = 0;
+elseif islogical(DEM.Z)
+ DEM.Z(~C.Z) = false;
+end
+
+

@GRIDobj/private/builtincaller.m

@@ -26,7 +26,7 @@
  validatealignment(varargin{isGRIDobj},varargin{~isGRIDobj});
  end
 
- if ref == 1;
+ if ref == 1
  OUT.Z = builtin(funname,varargin{isGRIDobj}.Z,varargin{~isGRIDobj});
  else
  OUT.Z = builtin(funname,varargin{~isGRIDobj},varargin{isGRIDobj}.Z);

@STREAMobj/STREAMobj.m

@@ -182,6 +182,7 @@
 
  function distance = get.distance(S)
  % [dynamic property] distance from outlet
+
  distance = zeros(numel(S.x),1);
  for r = numel(S.ix):-1:1
  distance(S.ix(r)) = distance(S.ixc(r)) + ...

@STREAMobj/maplateral.m

@@ -45,13 +45,15 @@
 % S = STREAMobj(FD,'minarea',1000);
 % S = klargestconncomps(S);
 % S = trunk(S);
-% g = maplateral(S,gradient8(DEM),100,{@min @max});
+% g = maplateral(S,gradient8(DEM),100,{@min @max @mean});
 % plotdz(S,DEM)
 % yyaxis right
-% plotdzshaded(S,g)
+% plotdzshaded(S,g(:,[1 2])
+% hold on
+% plotdz(S,g(:,3))
 % ylabel('Hillslope gradient [-]')
 %
-% See also: STREAMobj, SWATHobj
+% See also: STREAMobj, SWATHobj, STREAMobj/smooth
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
 % Date: 12. June, 2018

@STREAMobj/orientation.m

@@ -24,6 +24,7 @@
 % from top, wheras radians will be measured anti-clockwise from
 % the x-axis. 'cart' returns the cartesian coordinates of the
 % directional vectors.
+% 
 %
 % Output arguments
 %
@@ -71,5 +72,5 @@
 switch unit
  case 'degrees'
  theta = rad2deg(theta);
- theta = mod(-90-theta,360);
+ theta = mod(-270-theta,360);
 end

compilemexfiles.m

@@ -47,6 +47,7 @@
  for r = 1:numel(files)
  funname = files(r).name;
  mex('-largeArrayDims',funname)
+ % mex('-R2018a',funname)
  end
 catch err
  cd(oldFolder);

readme.md

@@ -93,6 +93,7 @@ resource for code and examples is the [TopoToolbox blog](https://topotoolbox.word
 - Documentation in the documentation browser
 - new function: ttcmap for creating nice colormaps for DEMs, particularly
  if DEMs include topography and bathymetry
+- FLOWobj2gradient renamed to gradient
 - new function: ttscm for access to scientific colormaps;
  see [Fabio Crameri's website](https://www.fabiocrameri.ch/colourmaps.php)
 - new function: mappingapp (lightweighed GUI for mapping points simultaneously in
@@ -106,6 +107,7 @@ resource for code and examples is the [TopoToolbox blog](https://topotoolbox.word
 - new function: STREAMobj/sinuosity
 - new function: STREAMobj/clean
 - new function: STREAMobj/nal2nal
+- new function: GRIDobj/clip
 - new function: GRIDobj/GRIDobj2im
 - new function: GRIDobj/getextent
 - new function: IOtools/readexample