Merge pull request #4 from wschwanghart/master

TopoToolbox 2.2

.gitignore

@@ -0,0 +1 @@
+*.asv

@FLOWobj/FLOWobj2GRIDobj.m

@@ -1,10 +1,10 @@
 function G = FLOWobj2GRIDobj(F)
 
-% create ESRI ArcGIS flow direction grid from FLOWobj
+%FLOWOBJ2GRIDOBJ create ESRI ArcGIS flow direction grid from FLOWobj
 %
 % Syntax
 %
-% G = FLOWobj2GRIDobj(F)
+% AF = FLOWobj2GRIDobj(FD)
 %
 % Description
 %
@@ -15,20 +15,30 @@
 % 
 % Input arguments
 %
-% F  instance of FLOWobj
+% FD instance of FLOWobj
 %
 % Output arguments
 %
-% G  instance of GRIDobj that stores the flow directions in the
+% AF instance of GRIDobj that stores the flow directions in the
 % format used by ESRI ArcGIS.
 %
+% Example
+%
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% FD = FLOWobj(DEM,'preprocess','c');
+% ArcFlow = FLOWobj2GRIDobj(FD);
+% imagesc(ArcFlow)
+%
 % 
 % See also: GRIDobj, GRIDobj2geotiff
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
 % Date: 5. October, 2013
 
-
+switch lower(F.type)
+ case {'multi' 'dinf'}
+ error('not possible for multiple flow directions')
+end
 
 G = copy2GRIDobj(F);
 

@FLOWobj/FLOWobj2M.m

@@ -1,6 +1,6 @@
 function M = FLOWobj2M(FD)
 
-% convert instance of FLOWobj to flow direction matrix 
+%FLOWOBJ2M convert instance of FLOWobj to flow direction matrix 
 %
 % Syntax
 %
@@ -11,9 +11,26 @@
 % FLOWobj2M converts an instance of FLOWobj to the flow direction
 % matrix M as used in TopoToolbox 1.
 %
+% Input arguments 
+%
+% FD FLOWobj
+%
+% Output arguments
+%
+% M sparse transfer matrix
+%
+% Example
+%
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% FD = FLOWobj(DEM,'preprocess','c');
+% M = FLOWobj2M(FD);
+% [x,y] = getcoordinates(DEM);
+% [x,y] = meshgrid(x,y);
+% gplot(M,[x(:) y(:)]) 
+%
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
-% Date: 23. February, 2013
+% Date: 18. August, 2017
 
 
 nrc = prod(FD.size);

@FLOWobj/FLOWobj2cell.m

@@ -1,6 +1,6 @@
 function [CFD,D,A,cfix] = FLOWobj2cell(FD,IX)
 
-% Return cell array of FLOWobjs for individual drainage basins
+%FLOWOBJ2CELL return cell array of FLOWobjs for individual drainage basins
 %
 % Syntax
 %
@@ -31,6 +31,13 @@
 %
 % Example
 %
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% FD = FLOWobj(DEM,'preprocess','c');
+% [CFD,~,a] = FLOWobj2cell(FD);
+% [~,ix] = max(a);
+% A = flowacc(CFD{ix});
+% imageschs(DEM,log(A))
+%
 % See also: FLOWobj, STREAMobj/STREAMobj2cell
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
@@ -49,6 +56,8 @@
  error('a forth output argument is only allowed for two input arguments');
 end
 
+
+
 FDcopy = FD;
 FDcopy.ix = [];
 FDcopy.ixc = [];

@FLOWobj/FLOWobj2gradient.m

@@ -1,6 +1,6 @@
 function G = FLOWobj2gradient(FD,DEM)
 
-% gradient along flow direction
+%FLOWOBJ2GRADIENT gradient along flow direction
 %
 % Syntax 
 %
@@ -13,6 +13,24 @@
 % when deriving the FLOWobj. FLOWobj2gradient returns an instance of
 % GRIDobj that contains the gradient along the flow paths.
 %
+% If FLOWobj has been derived using a multiple or Dinf flow direction
+% algorithm, then G is calculated as the weighted mean gradient.
+%
+% Input arguments
+%
+% FD FLOWobj
+% DEM digital elevation model (GRIDobj)
+%
+% Output arguments
+%
+% G along-flow gradient
+%
+% Example
+%
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% FD = FLOWobj(DEM,'preprocess','c');
+% G = FLOWobj2gradient(FD,DEM);
+% imageschs(DEM,G)
 %
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
@@ -21,5 +39,13 @@
 validatealignment(FD,DEM)
 d = getdistance(FD.ix,FD.ixc,FD.size,FD.cellsize);
 G = DEM;
-G.Z = zeros(FD.size,class(DEM.Z));
-G.Z(FD.ix) = (DEM.Z(FD.ix)-DEM.Z(FD.ixc))./d;
+switch FD.type
+ case 'single'
+ G.Z = zeros(FD.size,class(DEM.Z));
+ G.Z(FD.ix) = (DEM.Z(FD.ix)-DEM.Z(FD.ixc))./d;
+ otherwise
+ d = FD.fraction .* (DEM.Z(FD.ix) - DEM.Z(FD.ixc))./d;
+ g = accumarray(FD.ix,d,[prod(DEM.size) 1],@sum,zeros(1,'like',DEM.Z),false);
+ G.Z = reshape(g,DEM.size);
+end
+

@FLOWobj/dbentropy.m

@@ -0,0 +1,208 @@
+function [OUT,Pc] = dbentropy(FD,ix)
+
+%DBENTROPY entropy of drainage basin delineation
+%
+% Syntax
+%
+% [H,S] = dbentropy(M,siz,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.
+%
+% Input parameters
+%
+% FD multiple flow direction FLOWobj
+% IX linear indices into the DEM from which FD was derived.
+%
+% Output arguments
+% 
+% H Shannon Entropy grid (GRIDobj)
+% S structure array with numel(ix) elements with additional 
+% information to each outlet. S.P contains the probability grid
+% for each outlet. S.IX is the outlet's linear index. CA01-99 are
+% the error bounds for the drainage area.
+%
+% Example
+%
+% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+% DEM = fillsinks(DEM);
+% FD = FLOWobj(DEM,'multi');
+% IX = [107807 769639];
+% [H,SB] = dbentropy(FD,IX);
+% imageschs(DEM,H,'colormap',flipud(hot))
+%
+% Reference
+%
+% Schwanghart, W., Heckmann, T. (2012): Fuzzy delineation of drainage 
+% basins through probabilistic interpretation of diverging flow algorithms. 
+% Environmental Modelling and Software, 33, 106-113. 
+% [DOI: 10.1016/j.envsoft.2012.01.016]
+%
+% See also: drainagebasins, FLOWobj
+%
+% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
+% Date: 07. October, 2016
+
+M = FLOWobj2M(FD);
+siz = FD.size;
+
+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;
+ chunks = true;
+ else
+ chunks = false;
+ end
+
+ % In addition, suggest to output only some of the larger fuzzy drainage basins 
+ 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');
+ ix = ix(ixsort);
+
+ end
+ nrc = prod(siz);
+
+else
+ nrc = prod(siz);
+ nrix = numel(ix);
+
+ % create sinks at provided indices
+ SINKS = ones(nrc,1);
+ SINKS(nrix) = 0;
+ M = spdiags(SINKS,0,nrc,nrc)*M;
+
+ chunks = false;
+end
+
+
+
+% S is a matrix where each outlet has its own column. 
+S = sparse(ix,1:nrix,1,nrc,nrix);
+
+% B is the membership (partition) matrix (?) (fuzzy partition matrix). 
+% B contains the membership grades. (calculated by membership function below)
+
+if nargout >= 2;
+ Pc = struct('P',{},'IX',{},'CA50',{},'CA99',{},...
+ 'CA95',{},'CA66',{},'CA33',{},...
+ 'CA05',{},'CA01',{});
+ % linear index
+ pl = (1:prod(siz))';
+ % percentiles
+ percentiles = [.99 .95 .66 .50 .33 .05 .01];
+end
+
+
+if ~chunks
+ % solve system of equations to get probabilities
+ P = (speye(nrc)-M)\S;
+
+ if nargout >= 2
+
+ for r = 1:nrix;
+ Pc(r).P = full(reshape(P(:,r),siz));
+ Pc(r).IX = ix(r);
+ ptemp = sort(full(P(:,r)),'descend');
+ [~,m] = unique(ptemp);
+ perc = interp1(ptemp(m),pl(m),percentiles);
+
+ Pc(r).CA50 = perc(4);
+ Pc(r).CA99 = perc(1);
+ Pc(r).CA95 = perc(2);
+ Pc(r).CA66 = perc(3);
+ Pc(r).CA33 = perc(5);
+ Pc(r).CA05 = perc(6);
+ Pc(r).CA01 = perc(7);
+
+ Pc(r).percentiles = percentiles;
+ Pc(r).area = perc;
+
+
+ end
+ end
+
+ sP = sum(P,2);
+ sP = sparse(max(1-sP,0));
+
+ H = reshape(-full(sum(spfun(@(x) x.*log2(x),[P sP]),2)),siz);
+
+
+else
+ H = zeros(siz);
+ counter = 1;
+ sP = zeros(prod(siz),1);
+
+ for r = 1:ceil(nrix/chunksize);
+ cols = ((r-1)*chunksize + 1):r*chunksize;
+ cols(cols>nrix) = [];
+
+ P = (speye(nrc)-M)\S(:,cols);
+ sP = sum(P,2)+sP;
+ H = H + reshape(-full(sum(spfun(@(x) x.*log2(x),P),2)),siz);
+
+ if nargout >= 2
+
+ for rr = cols(1):cols(end);
+
+ if rr<=pout
+ Pc(rr).P = full(reshape(P(:,counter),siz));
+ ptemp = sort(Pc(rr).P(:),'descend');
+ else
+ Pc(rr).P = [];
+ ptemp = sort(full(P(:,counter)),'descend');
+ end
+
+ Pc(rr).IX = ix(rr);
+ [~,m] = unique(ptemp);
+ perc = interp1(ptemp(m),pl(m),percentiles);
+
+ Pc(rr).CA50 = perc(4);
+ Pc(rr).CA99 = perc(1);
+ Pc(rr).CA95 = perc(2);
+ Pc(rr).CA66 = perc(3);
+ Pc(rr).CA33 = perc(5);
+ Pc(rr).CA05 = perc(6);
+ Pc(rr).CA01 = perc(7);
+
+
+ % reset counter if larger than chunksize
+ if counter >= chunksize;
+ counter = 1;
+ else
+ counter = counter+1;
+ end
+
+
+ end
+ end
+ end 
+% H = H + reshape(1 - sP.*log2(sP),siz);
+end
+
+H = max(min(H,1),0);
+OUT = GRIDobj(FD);
+OUT.Z = H;
+OUT.name = 'fuzzy drainage basins';
+
+