few updates (e.g. fastscape)

@GRIDobj/GRIDobj2polygon.m

@@ -67,7 +67,7 @@
 % See also: bwboundaries, bwtraceboundary, regionprops
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
-% Date: 17. August, 2017
+% Date: 14. June, 2022
 
 
 
@@ -96,7 +96,7 @@
 waitb = p.Results.waitbar;
 
 % check underlying class of the grid
-if isfloat(DB.Z);
+if isfloat(DB.Z)
  writevalue = true;
  DB2 = GRIDobj(DB,'uint32');
  I = ~(isnan(DB.Z) | DB.Z == 0);
@@ -107,7 +107,11 @@
 end
 
 % identify regions and number of regions
-STATS = regionprops(uint32(DB.Z),'Area','PixelIdxList');
+if islogical(DB.Z)
+ STATS = regionprops(DB.Z,'Area','PixelIdxList');
+else
+ STATS = regionprops(uint32(DB.Z),'Area','PixelIdxList');
+end
 ndb = numel(STATS); 
 
 % go through all regions

@PPS/extractvaluesaroundpoints.m

@@ -1,12 +1,12 @@
-function [v,nIX] = extractvaluesaroundpoints(P,z,varargin)
+function [v,nIX,CP] = extractvaluesaroundpoints(P,z,varargin)
 
 %EXTRACTVALUESAROUNDPOINTS Extract values around points 
 %
 % Syntax
 %
 % v = extractvaluesaroundpoints(P,z)
 % v = extractvaluesaroundpoints(p,z,'pn',pv,...)
-% [v,nIX] = ...
+% [v,nIX,CP] = ...
 %
 % Description
 %
@@ -26,9 +26,18 @@
 % 
 % Parameter name/value pairs
 %
-% 'direction' {'both'}, 'upstream' or 'downstream'.
-% 'dfrompoint' distance from points
-% 'aggfun' anonymous aggregation function. {@mean}
+% 'direction' {'both'}, 'upstream', 'downstream' or 'bothundirected'
+% 'both' extracts the stream network in upstream and
+% downstream direction. It won't include downstream
+% tributaries, however. 'bothundirected' will include all
+% stream segments, irrespective of their direction.
+% 'dfrompoint' distance from points. If 'direction' is 'both', then
+% 'dfrompoint' accepts a two-element vector where there
+% first element refers to the distance downstream from
+% the point, and the second element refers to the
+% distance upstream from the point.
+% 'aggfun' anonymous aggregation function which takes a vector and 
+% returns a scalar {@mean}. 
 % 'distance' By default, distances between are calculated as the
 % euclidean distance between stream nodes. Yet, other
 % distance metrics (e.g. chi) can be used as well. Note
@@ -38,8 +47,9 @@
 % Output arguments
 %
 % v extracted value
-% nIX cell array of linear indices into node-attribute list for all 
-% neighbors of all points in P.
+% nIX cell array of linear indices into the grid.
+% CP cell array of PPS objects, one for each point, with the river
+% network extracted. 
 % 
 % Example
 %
@@ -51,7 +61,7 @@
 p = inputParser;
 p.FunctionName = 'extractvaluesaroundpoints';
 
-addParameter(p,'direction','both',@(x) ischar(validatestring(x,{'both','upstream','downstream'})))
+addParameter(p,'direction','both',@(x) ischar(validatestring(x,{'both','upstream','downstream','bothundirected'})))
 addParameter(p,'dfrompoint',10*P.S.cellsize)
 addParameter(p,'distance',P.S.distance)
 addParameter(p,'aggfun',@mean)
@@ -79,19 +89,72 @@
  G = digraph(P.S.ix, P.S.ixc,w);
  case 'upstream'
  G = digraph(P.S.ixc, P.S.ix, w);
- case 'both'
+ case 'bothundirected'
  G = graph([P.S.ix; P.S.ixc],[P.S.ixc; P.S.ix], [w;w]);
+ case 'both'
+ % downstream
+ G1 = digraph(P.S.ix, P.S.ixc,w);
+ d1 = p.Results.dfrompoint(1); 
+ % upstream
+ G2 = digraph(P.S.ixc, P.S.ix, w);
+ if numel(p.Results.dfrompoint) == 2
+ d2 = p.Results.dfrompoint(2);
+ else
+ d2 = d1;
+ end
 end
 
 % Calculate nearest neighbor indices using cellfun
 pIX = num2cell(P.PP);
 
-if nargout == 2
- nIX = cellfun(@(s) nearest(G,s,p.Results.dfrompoint),pIX,'UniformOutput',false);
+if nargout >= 2
+ switch direction 
+ case 'both'
+
+ nIX = cellfun(@(s) ...
+ [nearest(G1,s,d1);...
+ s; ...
+ nearest(G2,s,d2)],...
+ pIX,'UniformOutput',false);
+
+
+ otherwise
+ nIX = cellfun(@(s) [nearest(G,s,p.Results.dfrompoint(1));s],pIX,'UniformOutput',false);
+ end
  v = cellfun(@(ix) p.Results.aggfun(z(ix)),nIX);
+
 elseif nargout == 1
- v = cellfun(@(s) p.Results.aggfun(z(nearest(G,s,p.Results.dfrompoint))),...
+ switch direction 
+ case 'both'
+ v = cellfun(@(s) p.Results.aggfun(z(...
+ [nearest(G1,s,d1);...
+ s; ...
+ nearest(G2,s,d2)])),...
  pIX,'UniformOutput',true);
+
+ otherwise
+ v = cellfun(@(s) p.Results.aggfun(z([nearest(G,s,p.Results.dfrompoint(1));s])),...
+ pIX,'UniformOutput',true);
+ end
 end
 
-
+% Extract stream network around points
+CP = cell(npoints(P),1);
+if nargout >= 3
+ for r = 1:npoints(P)
+ s = logical(getnal(P.S));
+ s(nIX{r}) = true;
+
+ if isempty(P.z)
+ S = subgraph(P.S,s);
+ CP{r} = PPS(S,'PP',P.S.IXgrid(P.PP(r)));
+ else
+ [S,locb] = subgraph(P.S,s);
+ CP{r} = PPS(S,'PP',P.S.IXgrid(P.PP(r)),'z',P.z(locb));
+ end
+ end
+end
+
+if nargout >= 2
+ nIX = cellfun(@(ix) P.S.IXgrid(ix),nIX,'UniformOutput',false);
+end

@STREAMobj/fastscape.m

@@ -45,6 +45,12 @@
 % dt time step length [y] (default is 1000)
 % plot true or false (default is true)
 % ploteach plot each time step (default is 10)
+% plotchi {false} or true. If true, the plot will show the
+% chi-transformed profile.
+% gifname By default empty, but if filename is provided, a gif file
+% will be written to the disk. 
+% gifopts Structure array with gif options for the function gif
+% (see help gif)
 %
 % Example
 %
@@ -101,13 +107,23 @@
 % Hergarten, S., 2002. Self organised criticality in Earth Systems. 
 % Springer, Heidelberg.
 %
-% See also: STREAMobj
+% Campforts, B., Schwanghart, W., Govers, G. (2017): Accurate 
+% simulation of transient landscape evolution by eliminating numerical 
+% diffusion: the TTLEM 1.0 model. Earth Surface Dynamics, 5, 47-66. 
+% [DOI: 10.5194/esurf-5-47-2017]
+%
+% See also: STREAMobj, gif
 %
 % Authors: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de) and
 % Benjamin Campforts.
 % Date: 28. April, 2022
 
 
+gifopts.DelayTime = 1/15;
+gifopts.LoopCount = inf;
+gifopts.frame = gcf;
+gifopts.overwrite = false;
+
 %% Input parsing
 p = inputParser;
 addRequired(p,'S',@(x) isa(x,'STREAMobj'))
@@ -125,6 +141,9 @@
 addParameter(p,'bctype','rate')
 addParameter(p,'plot',true)
 addParameter(p,'ploteach',10)
+addParameter(p,'plotchi',0)
+addParameter(p,'gifname','')
+addParameter(p,'gifopts',gifopts)
 parse(p,S,z,a,varargin{:})
 
 % STREAMobj
@@ -150,6 +169,11 @@
 % Plot?
 plotit = p.Results.plot;
 ploteach = p.Results.ploteach;
+plotchi = p.Results.plotchi;
+% Write gif
+writegif = p.Results.gifname;
+gifopts = p.Results.gifopts;
+
 
 %FASTSCAPE1D 1D implementation of Braun and Willett 2013 implicit scheme
 ix = S.ix;
@@ -158,6 +182,7 @@
 dx_ixcix = d(ix)-d(ixc);
 
 % Calculate K*A^m
+ar = a;
 a = k.*(a.^m);
 
 % calculate timesteps
@@ -202,8 +227,20 @@
 end
 
 if plotit
- hh = plotdz(S,z,'color','r');
+ if plotchi == 0
+ d = S.distance;
+ elseif plotchi == 1
+ d = chitransform(S,ar,"mn",m/n);
+ end
+ hh = plotdz(S,z,'color','r','distance',d); 
  hold on
+ if plotchi
+ xlabel('\chi [m]')
+ end
+
+ if ~isempty(writegif)
+ gif(writegif,gifopts)
+ end
 end
 
 t = 0;
@@ -231,7 +268,10 @@
  delete(h(1))
  h(1:9) = h(2:10);
  end 
- h(plotcounter) = plotdz(S,z,'color','k');
+ h(plotcounter) = plotdz(S,z,'color','k','distance',d);
+ if plotchi
+ xlabel('\chi [m]')
+ end
  if plotcounter > 2
  for rr = 1:plotcounter
  h(rr).Color = [0 0 0 rr/(plotcounter*2)];
@@ -244,6 +284,10 @@
  end
 
  drawnow
+ if ~isempty(writegif)
+ gif
+ end
+
  end
  end
 
@@ -287,7 +331,12 @@
  delete(h(1))
  h(1:9) = h(2:10);
  end 
- h(plotcounter) = plotdz(S,z,'color','k');
+ h(plotcounter) = plotdz(S,z,'color','k','distance',d);
+ if plotchi
+ xlabel('\chi [m]')
+ end
+
+
  if plotcounter > 2
  for rr = 1:plotcounter
  h(rr).Color = [0 0 0 rr/(plotcounter*2)];
@@ -300,6 +349,9 @@
  end
 
  drawnow
+ if ~isempty(writegif)
+ gif
+ end
  end
  end
 end

@STREAMobj/modify.m

@@ -173,7 +173,7 @@
 addParamValue(p,'rmconncomps_ch',[],@(x) isnumeric(x) && x>=0 && isscalar(x));
 addParamValue(p,'rmupstreamtoch',[],@(x) isa(x,'STREAMobj'));
 addParamValue(p,'fromch2IX',[]);
-addParamValue(p,'rmnodes',[],@(x) isa(x,'STREAMobj'));
+addParamValue(p,'rmnodes',[]);
 addParamValue(p,'clip',[],@(x) (isnumeric(x) && size(x,2)==2 && size(x,1)>2) || isa(x,'GRIDobj'));
 addParamValue(p,'nal',[],@(x) isnal(S,x));
 
@@ -294,7 +294,7 @@
 % end
  end
 
- I = false(size(S.x));
+
  I = II;
  for r = 1:numel(S.ix)
  I(S.ixc(r)) = II(S.ix(r)) || I(S.ix(r)) || I(S.ixc(r));
@@ -378,8 +378,13 @@
  I = ismember(cc,md);
 
 elseif ~isempty(p.Results.rmnodes)
- I = ~ismember(S.IXgrid,p.Results.rmnodes.IXgrid);
-
+ % check if empty S is supplied
+ if isempty(p.Results.rmnodes.IXgrid)
+ I = getnal(S) == 0;
+ else
+ I = ~ismember(S.IXgrid,p.Results.rmnodes.IXgrid);
+ end
+
 elseif ~isempty(p.Results.nal)
  I = p.Results.nal;
 
@@ -397,7 +402,7 @@
 
 elseif ~isempty(p.Results.interactive)
 %% interactive 
- figure
+%  figure
  plot(S,'k'); axis equal
  ax = gca;
  % expand axes
@@ -461,14 +466,25 @@
  switch meth
  case 'polyselect'
  title('create a polygon and double-click to finalize')
+ try
+ hp = drawpolygon;
+ pos = hp.Position;
+ catch
  hp = impoly;
  pos = wait(hp);
+ end
  pos(end+1,:) = pos(1,:);
  case 'ellipseselect'
  title('create a ellipse and double-click to finalize')
+ try
+ hp = drawellipse;
+ pos = hp.Vertices;
+ catch
  hp = imellipse;
- wait(hp);
+ pos = wait(hp);
  pos = getVertices(hp);
+ end
+
  pos(end+1,:) = pos(1,:);
  case 'rectselect'
  title('create a rectangle and double-click to finalize')