update

@GRIDobj/cellarea.m

@@ -1,6 +1,6 @@
 function CA = cellarea(DEM,unit)
 
-% calculate cell areas of a GRIDobj in geographic coordinate system
+%CELLAREA calculate cell areas of a GRIDobj in geographic coordinate system
 %
 % Syntax
 %
@@ -26,16 +26,20 @@
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
 % Date: 20. May, 2016
 
-if isempty(DEM.georef)
- error('no coordinate system defined')
+
+try 
+ GST = DEM.georef.SpatialRef.CoordinateSystemType;
+catch
+ GST = 'geographic';
 end
 
-switch DEM.georef.SpatialRef.CoordinateSystemType
+switch GST
  case 'geographic'
  otherwise
  error('DEM has a projected coordinate system')
 end
 
+
 switch lower(unit)
  case 'm'
  scale = 1e6;

@GRIDobj/project.m

@@ -0,0 +1,196 @@
+function DEMr = project(SOURCE,TARGET,varargin)
+
+%PROJECT transforms a GRIDobj between projected coordinate systems
+%
+% Syntax
+%
+% OUT = project(SOURCE,TARGET)
+% OUT = project(SOURCE,TARGET,pn,pv,...)
+%
+% Description
+%
+% project reprojects the GRIDobj 'SOURCE' to have the same projection
+% as the GRIDobj 'TARGET'. Both GRIDobj's need to have an mstruct in
+% their 'georef' fields. If SOURCE does not, then it is assumed to have
+% a geographic coordinate system (horizontal WGS84 datum). The function
+% transforms the projected coordinates first to geographic coordinates
+% and then back to projected coordinates using the functions mfwdtran
+% and minvtran (both part of the Mapping Toolbox). The spatial
+% resolution of the output (OUT) will be the same as of TARGET, unless
+% set differently by the optional variable res.
+%
+% Input arguments
+%
+% SOURCE instance of GRIDobj that shall be transformed
+% TARGET instance of GRIDobj with the target projection
+%
+% Parameter name/value pairs
+% 
+% 'res' scalar
+% spatial resolution. Default is TARGET.cellsize
+% 'method' string
+% 'bilinear' (default),'bicubic' or 'nearest'
+% 'fillvalue' scalar
+% nan (for single and double grids). Otherwise 0.
+% 'align' logical scalar
+% true (default) or false. If true, the OUT grid will be
+% spatially aligned with TARGET. This means, they have
+% the same extent and spatial alignment of cells. If
+% 'align, true then setting the resolution 'res' will
+% have no effect.
+% 
+% Output arguments
+%
+% OUT instance of GRIDobj
+%
+%
+% See also: GRIDobj, reproject2utm
+%
+% Author: Dirk Scherler (scherler[at]gfz-potsdam.de) and Wolfgang
+% Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
+% Date: 2. December, 2018
+
+
+p = inputParser;
+p.FunctionName = 'GRIDobj/project';
+addParamValue(p,'res', TARGET.cellsize,@(x) isscalar(x));
+addParamValue(p,'align',true,@(x) isscalar(x));
+addParamValue(p,'method','bilinear',@(x) ischar(x));
+addParamValue(p,'fillvalue',nan,@(x) isscalar(x));
+parse(p,varargin{:});
+
+% Get mstruct of SOURCE
+try
+ mstruct_s = SOURCE.georef.mstruct;
+ %cstype_s = SOURCEDEM.georef.SpatialRef.CoordinateSystemType;
+ pr2pr = true;
+catch
+ % there is either no mstruct available or the coordinate system is
+ % unknown. We assume that the source coordinate system is geographic
+ % (WGS84 Datum)
+ pr2pr = false;
+end
+
+% Get mstruct of TARGET
+mstruct_t = TARGET.georef.mstruct;
+%cstype_t = TARGETDEM.georef.SpatialRef.CoordinateSystemType;
+
+if ~p.Results.align 
+
+ [x0,y0] = getcoordinates(SOURCE);
+ x0 = x0(:);
+ y0 = y0(:);
+
+ % Calculate bounds of the reprojected DEM
+ x = [x0; x0; repmat(x0(1),numel(y0),1); repmat(x0(end),numel(y0),1)];
+ y = [repmat(y0(1),numel(x0),1); repmat(y0(end),numel(x0),1); y0; y0]; 
+
+ if pr2pr
+ [lat,lon] = minvtran(mstruct_s,x,y);
+ [x,y] = mfwdtran(mstruct_t,lat,lon);
+ else
+ [x,y] = mfdwtran(mstruct_t,y,x);
+ end
+
+ tlims([1 2]) = [min(x),max(x)];
+ tlims([3 4]) = [min(y),max(y)];
+
+ res = p.Results.res;
+else
+
+ [x0,y0] = getcoordinates(SOURCE);
+ [x,y] = getcoordinates(TARGET);
+ tlims([1 2]) = [min(x), max(x)];
+ tlims([3 4]) = [min(y),max(y)];
+ res = TARGET.cellsize;
+
+end
+
+% Limits of source grid
+slims([1 2]) = [min(x0) max(x0)];
+slims([3 4]) = [min(y0) max(y0)];
+
+% get fillvalue
+fillval = p.Results.fillvalue;
+if isinteger(SOURCE.Z)
+ fillval = zeros(1,class(SOURCE.Z));
+elseif islogical(SOURCE.Z)
+ fillval = false;
+else
+ fillval = cast(fillval,class(SOURCE.Z));
+end
+fillval = double(fillval);
+
+% check method
+meth = validatestring(p.Results.method,{'bilinear','linear','nearest','bicubic'});
+
+% Prepare tform for the image transform
+if pr2pr
+ T = maketform('custom', 2, 2, @FWDTRANSpr2pr, @INVTRANSpr2pr, []);
+else
+ T = maketform('custom', 2, 2, @FWDTRANSgcs2pr, @INVTRANSgcs2pr, []);
+end
+
+% Calculate image transform
+[Znew,xdata,ydata] = imtransform(flipud(SOURCE.Z),T,meth,...
+ 'Xdata',tlims([1 2]),'Ydata',tlims([3 4]),...
+ 'Udata',slims([1 2]),'Vdata',slims([3 4]),...
+ 'XYScale',[res res],'Fillvalues',fillval);
+
+Znew = flipud(Znew); 
+
+if p.Results.align
+ % If aligned, this is easy. The transformed GRIDobj will be perfectly
+ % aligned with TARGET
+ DEMr = TARGET;
+ DEMr.Z = Znew;
+else
+ % We have calculated the imtransform with 'ColumnsStartFrom' south. 
+ % GRIDobjs use 'ColumnsStartFrom' north
+
+ xnew = cumsum([xdata(1) repmat(res,1,size(Znew,2)-1)]);
+ ynew = flipud(cumsum([ydata(1) repmat(res,1,size(Znew,1)-1)])');
+
+ % Construct GRIDobj
+ DEMr = GRIDobj(xnew,ynew,Znew);
+ % figure, subplot(1,2,1), imagesc(DEMr), subplot(1,2,2), imagesc(TARGETDEM)
+
+ % Include geospatial and other information
+ R = refmatToMapRasterReference(DEMr.refmat,DEMr.size);
+ DEMr.georef.SpatialRef = R;
+ DEMr.georef.mstruct = TARGET.georef.mstruct;
+ DEMr.georef.GeoKeyDirectoryTag = TARGET.georef.GeoKeyDirectoryTag;
+ DEMr.zunit = SOURCE.zunit;
+ DEMr.xyunit = SOURCE.xyunit;
+end
+
+DEMr.name = [SOURCE.name ' (projected)'];
+
+
+% Transformation functions for imtransform (projected --> projected)
+ function x = FWDTRANSpr2pr(u,~)
+ [tlat,tlon] = minvtran(mstruct_s,u(:,1),u(:,2));
+ [tx,ty] = mfwdtran(mstruct_t,tlat,tlon);
+ x = [tx ty]; 
+ end
+
+
+ function u = INVTRANSpr2pr(x,~)
+ [tlat,tlon] = minvtran(mstruct_t,x(:,1),x(:,2));
+ [tx,ty] = mfwdtran(mstruct_s,tlat,tlon);
+ u = [tx ty];
+ end
+
+% Transformation functions for imtransform (gcs --> projected)
+ function x = FWDTRANSgcs2pr(u,~)
+ [tx,ty] = mfwdtran(mstruct_t,u(:,2),u(:,1));
+ x = [tx ty]; 
+ end
+
+
+ function u = INVTRANSgcs2pr(x,~)
+ [tlat,tlon] = minvtran(mstruct_t,x(:,1),x(:,2));
+ u = [tlon tlat];
+ end
+
+end

@STREAMobj/labelreach.m

@@ -49,37 +49,82 @@
 p.FunctionName = 'STREAMobj/labelreach';
 addParamValue(p,'seglength',inf,@(x) isscalar(x) && x>0);
 addParamValue(p,'shuffle',false,@(x) isscalar(x));
+addParamValue(p,'exact',false,@(x) isscalar(x));
 parse(p,varargin{:});
 
-label = streampoi(S,{'bconfl','out'},'logical');
-label = cumsum(label).*label;
 
-ix = S.ix;
-ixc = S.ixc;
-for r = numel(ix):-1:1
- if label(ix(r)) == 0
- label(ix(r)) = label(ixc(r));
+if ~p.Results.exact
+ % The first approach will label all nodes in the stream network.
+ % It will split reaches according to segment length but also at
+ % confluences. This will result in segment lengths that are not exactly
+ % the provided value seglength.
+ label = streampoi(S,{'bconfl','out'},'logical');
+ label = cumsum(label).*label;
+
+ ix = S.ix;
+ ixc = S.ixc;
+ for r = numel(ix):-1:1
+ if label(ix(r)) == 0
+ label(ix(r)) = label(ixc(r));
+ end
  end
-end
-
-% there may be nodes that are neither receivers or givers. These nodes
-% still have the label zero. 
-maxlabel = max(label);
-iszero = label==0;
-label(iszero) = (1:nnz(iszero))+maxlabel;
-
-
-if ~isinf(p.Results.seglength)
- maxdist = p.Results.seglength;
- cs = S.cellsize;
- dist = S.distance;
- label2 = accumarray(label,(1:numel(S.x))',[max(label) 1],@(x) {split(x)});
- labeltemp = zeros(size(label));
- for r = 1:max(label)
- labeltemp(label2{r}(:,2)) = label2{r}(:,1);
+
+ % there may be nodes that are neither receivers or givers. These nodes
+ % still have the label zero.
+ maxlabel = max(label);
+ iszero = label==0;
+ label(iszero) = (1:nnz(iszero))+maxlabel;
+
+
+ if ~isinf(p.Results.seglength)
+ maxdist = p.Results.seglength;
+ cs = S.cellsize;
+ dist = S.distance;
+ label2 = accumarray(label,(1:numel(S.x))',[max(label) 1],@(x) {split(x)});
+ labeltemp = zeros(size(label));
+ for r = 1:max(label)
+ labeltemp(label2{r}(:,2)) = label2{r}(:,1);
+ end
+ [~,~,label] = unique([label labeltemp],'rows');
+ end
+else
+ % If we choose the exact approach, then labelreach will not split at
+ % confluences. Parts close to channelheads that have length less than
+ % the value of seglength will get the label zero.
+ [CS,locS] = STREAMobj2cell(S,'tributaries');
+ labeltrib = zeros(size(S.x));
+ for r = 1:numel(locS)
+ labeltrib(locS{r}) = r;
  end
- [~,~,label] = unique([label labeltemp],'rows');
+
+ label = zeros(size(S.x));
+ Clabel = cell(size(CS));
+
+ for iter = 1:numel(CS)
+ d = CS{iter}.distance;
+ d = mod(d,p.Results.seglength);
+ I = double(gradient(CS{iter},d) < 0);
+ I(streampoi(CS{iter},'outlet','logical')) = 1;
+
+ ix = CS{iter}.ix;
+ ixc = CS{iter}.ixc;
+
+ for r = numel(ix):-1:1
+ I(ix(r)) = I(ixc(r)) + I(ix(r));
+ end
+ Clabel{iter} = I;
+ end
+
+ for r = numel(CS):-1:1
+ label(locS{r}) = Clabel{r};
+ end
+
+ % get unique labels
+ [~,~,label] = unique([label labeltrib],'rows');
+
 end
+
+
 
 %% Shuffle labels
 if p.Results.shuffle

@STREAMobj/split.m

@@ -6,7 +6,8 @@
 %
 % Ss = split(S)
 % Ss = split(S,ix)
-% Ss = split(S,ix,removeedge)
+% Ss = split(S,C)
+% Ss = split(S,...,removeedge)
 %
 % Description
 %
@@ -24,6 +25,8 @@
 % S STREAMobj
 % ix linear index into DEM. The indexed pixels must be located on 
 % the channel network, e.g. ismember(ix,S.IXgrid).
+% C GRIDobj. The GRIDobj should be a label grid. The stream network
+% is split where it traverses different regions. 
 % removeedge {'incoming'} or 'outgoing'. 'incoming' removes the edge or
 % edges between ix and its upstream neighbor(s). 'outgoing'
 % removes the downstream edge.
@@ -52,14 +55,19 @@
 % STREAMobj/STREAMobj2cell
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
-% Date: 15. May, 2017
+% Date: 3. December, 2018
 
 
 narginchk(1,3)
 if nargin == 1
  V = streampoi(S,'confluences','logical');
 else
- V = ismember(S.IXgrid,varargin{1});
+ if isa(varargin{1},'GRIDobj')
+ d = gradient(S,varargin{1});
+ V = d~=0;
+ else 
+ V = ismember(S.IXgrid,varargin{1});
+ end
 end
 
 if nargin <= 2
@@ -76,4 +84,6 @@
 end
 
 S.ix(I) = [];
-S.ixc(I) = [];
+S.ixc(I) = [];
+
+S = clean(S);