few updates

@FLOWobj/upslopestats.m

@@ -68,7 +68,8 @@
 if nargin == 2
  meth = 'mean';
 else
- meth = validatestring(meth,{'mean','std','var','min','max','sum'});
+ meth = validatestring(meth,{'mean','std','var','min','max','sum',...
+ 'nanmean'});
 end
 
 if nargin == 4
@@ -93,6 +94,15 @@
  % Sum
  S = flowacc(FD,VAR);
  S = S.Z;
+ case 'nanmean'
+ % Count
+ I = isnan(VAR);
+ n = flowacc(FD,+(~I));
+ n = n.Z;
+ % Sum
+ VAR(I) = 0;
+ S = flowacc(FD,VAR);
+ S = S.Z;
 end
 
 
@@ -105,6 +115,9 @@
  case 'mean'
  % Sum/Count
  VAR = S./n;
+ case 'nanmean'
+ VAR = S./n;
+ VAR(n==0) = nan; 
  case {'std','var'}
  % formula for calculating an unbiased estimate of the population variance
  % ! may be very much affected by round-off errors, arithmetic

@GRIDobj/filter.m

@@ -38,7 +38,9 @@
 
 
 
-validmethods = {'mean','average','median','sobel','scharr','wiener'};
+validmethods = {'mean','average','median',...
+ 'sobel','scharr','wiener',...
+ 'std'};
 
 % Parse inputs
 p = inputParser;
@@ -94,6 +96,8 @@
 
  case 'wiener'
  dem = wiener2(dem,ws);
+ case 'std'
+ dem = stdfilt(dem,true(p.Results.kernel(1),p.Results.kernel(2)));
 end
 
 % and set nans at their previous position

@GRIDobj/ne.m

@@ -0,0 +1,5 @@
+function OUT = ne(varargin)
+
+funname = 'ne';
+narginchk(2,2)
+OUT = builtincaller(funname,varargin{:});

@GRIDobj/project.m

@@ -22,12 +22,14 @@
 % Input arguments
 %
 % SOURCE instance of GRIDobj that shall be transformed
-% TARGET instance of GRIDobj with the target projection
+% TARGET instance of GRIDobj with the target projection or 
+% mapping projection structure (mstruct)
 %
 % Parameter name/value pairs
 % 
 % 'res' scalar
-% spatial resolution. Default is TARGET.cellsize
+% spatial resolution. Default is TARGET.cellsize. If
+% target is an mstruct, res must be supplied. 
 % 'method' string
 % 'bilinear' (default),'bicubic' or 'nearest'
 % 'fillvalue' scalar
@@ -51,9 +53,20 @@
 % Date: 2. December, 2018
 
 
+narginchk(2,inf)
+
+if isa(TARGET,'GRIDobj')
+ cs = TARGET.cellsize;
+ targetisGRIDobj = true;
+elseif isstruct(TARGET)
+ % target is a mapping structure
+ cs = [];
+ targetisGRIDobj = false;
+end
+
 p = inputParser;
 p.FunctionName = 'GRIDobj/project';
-addParamValue(p,'res', TARGET.cellsize,@(x) isscalar(x));
+addParamValue(p,'res', cs,@(x) isscalar(x));
 addParamValue(p,'align',true,@(x) isscalar(x));
 addParamValue(p,'method','bilinear',@(x) ischar(x));
 addParamValue(p,'fillvalue',nan,@(x) isscalar(x));
@@ -62,6 +75,9 @@
 % Get mstruct of SOURCE
 try
  mstruct_s = SOURCE.georef.mstruct;
+ if isempty(mstruct_s)
+ error('error')
+ end
  %cstype_s = SOURCEDEM.georef.SpatialRef.CoordinateSystemType;
  pr2pr = true;
 catch
@@ -72,10 +88,22 @@
 end
 
 % Get mstruct of TARGET
-mstruct_t = TARGET.georef.mstruct;
+if targetisGRIDobj
+ mstruct_t = TARGET.georef.mstruct;
+ cs = p.Results.res;
+else
+ mstruct_t = TARGET;
+ cs = p.Results.res;
+ if isempty(cs)
+ error('Spatial resolution of the output GRIDobj must be set.');
+ end
+end
+
+
+
 %cstype_t = TARGETDEM.georef.SpatialRef.CoordinateSystemType;
 
-if ~p.Results.align   
+if ~p.Results.align || ~targetisGRIDobj 
 
  [x0,y0] = getcoordinates(SOURCE);
  x0 = x0(:);
@@ -89,7 +117,7 @@
  [lat,lon] = minvtran(mstruct_s,x,y);
  [x,y] = mfwdtran(mstruct_t,lat,lon);
  else
- [x,y] = mfdwtran(mstruct_t,y,x);
+ [x,y] = mfwdtran(mstruct_t,y,x);
  end
 
  tlims([1 2]) = [min(x),max(x)];
@@ -139,7 +167,7 @@
 
 Znew = flipud(Znew); 
 
-if p.Results.align
+if p.Results.align && targetisGRIDobj
  % If aligned, this is easy. The transformed GRIDobj will be perfectly
  % aligned with TARGET
  DEMr = TARGET;
@@ -158,8 +186,13 @@
  % 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.georef.mstruct = mstruct_t;
+ if targetisGRIDobj
+ DEMr.georef.GeoKeyDirectoryTag = TARGET.georef.GeoKeyDirectoryTag;
+ else
+ DEMr.georef.GeoKeyDirectoryTag = [];
+ warning('Could not set GeoKeyDirectoryTag in output GRIDobj')
+ end
  DEMr.zunit = SOURCE.zunit;
  DEMr.xyunit = SOURCE.xyunit;
 end

@GRIDobj/surf.m

@@ -55,19 +55,45 @@
 % Date: 30. January, 2013
 
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% change threshold (maximum number of rows or cols) here 
-maxrowsorcols = 2000;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Parse inputs
+p = inputParser;
+p.KeepUnmatched = true;
+addOptional(p,'A',[]);
+addParameter(p,'maxrowsorcols',2000)
+addParameter(p,'exaggerate',1);
+addParameter(p,'block',false);
+addParameter(p,'baselevel',[]);
+addParameter(p,'sea',false);
+addParameter(p,'sealevel',0);
+addParameter(p,'seaalpha',0.5);
+
+% Parse
+parse(p,varargin{:});
+
+% create variables
+maxrowsorcols = p.Results.maxrowsorcols;
+exagg = p.Results.exaggerate;
+block = p.Results.block;
+baselevel = p.Results.baselevel;
+
+if block
+ minz = min(DEM);
+ maxz = max(DEM);
+ baselevel = minz-(maxz-minz)*0.2;
+end 
+sea = p.Results.sea;
+sealevel = p.Results.sealevel;
+seaalpha = p.Results.seaalpha;
+%%
 
 if max(DEM.size)>maxrowsorcols
 
  resamplecellsize = max(DEM.size)/(maxrowsorcols-1) * DEM.cellsize;
  DEM = resample(DEM,resamplecellsize);
 
- if ~isempty(varargin) && isa(varargin{1},'GRIDobj')
- A = varargin{1};
- varargin(1) = [];
+ if ~isempty(p.Results.A)
+ A = p.Results.A;
  A = resample(A,DEM);
  overlay = true;
  else
@@ -76,64 +102,30 @@
 
 
 else
- if ~isempty(varargin) && isa(varargin{1},'GRIDobj')
- A = varargin{1};
- varargin(1) = [];
+ if ~isempty(p.Results.A)
+ A = p.Results.A;
  validatealignment(DEM,A);
  overlay = true;
  else
  overlay = false;
  end
 end
 
-% go through varargin to find 'exaggerate'
-TF = strcmpi('exaggerate',varargin);
-ix = find(TF);
-if ~isempty(ix)
- exagg = varargin{ix+1};
- varargin([ix ix+1]) = [];
-else
- exagg = 1;
-end
-
-% go through varargin to find 'block'
-TF = strcmpi('block',varargin);
-ix = find(TF);
-if ~isempty(ix)
- block = varargin{ix+1};
- varargin([ix ix+1]) = [];
-else
- block = false;
-end
-
-% go through varargin to find 'baselevel'
-TF = strcmpi('baselevel',varargin);
-ix = find(TF);
-if ~isempty(ix)
- baselevel = varargin{ix+1};
- varargin([ix ix+1]) = [];
-else
- baselevel = min(DEM)*.8;
-end
 
-% go through varargin to find 'sea'
-TF = strcmpi('sea',varargin);
-ix = find(TF);
-if ~isempty(ix)
- sea = varargin{ix+1};
- varargin([ix ix+1]) = [];
-else
- sea = false;
-end
+[x,y] = refmat2XY(DEM.refmat,DEM.size);
 
+% Create cellarray from p.Unmatched
+pn = fieldnames(p.Unmatched);
+pv = struct2cell(p.Unmatched);
 
-
-[x,y] = refmat2XY(DEM.refmat,DEM.size); 
+pnpv = [pn pv];
+pnpv = pnpv';
+pnpv = pnpv(:)';
 
 if overlay
- h = surf(x,y,double(DEM.Z),double(A.Z),varargin{:});
+ h = surf(x,y,double(DEM.Z),double(A.Z),pnpv{:});
 else
- h = surf(x,y,double(DEM.Z),varargin{:});
+ h = surf(x,y,double(DEM.Z),pnpv{:});
 end
 
 exaggerate(gca,exagg);
@@ -194,55 +186,54 @@
  end
  facecolor = [1 1 1];
  facecolordark = [0.8 .8 .8];
- seaalpha = .5;
- sealevel = 0;
 
  xp = x(:);
  xp = [xp(1); xp; xp(end:-1:1)];
  yp = repmat(y(1),size(xp));
  zp = DEM.Z(1,:);
  zp = min(zp(:),0);
  zp = double([sealevel;zp;repmat(sealevel,size(zp))]);
- p(1) = patch(xp,yp,zp,facecolor,'EdgeColor','none','FaceColor',facecolor,'FaceAlpha',seaalpha);
+ pa(1) = patch(xp,yp,zp,facecolor,'EdgeColor','none','FaceColor',facecolor,'FaceAlpha',seaalpha);
 
  yp = repmat(y(end),size(xp));
  zp = DEM.Z(end,:);
  zp = min(zp(:),0);
  zp = double([sealevel;zp;repmat(sealevel,size(zp))]);
- p(2) = patch(xp,yp,zp,facecolor,'EdgeColor','none','FaceColor',facecolor,'FaceAlpha',seaalpha);
+ pa(2) = patch(xp,yp,zp,facecolor,'EdgeColor','none','FaceColor',facecolor,'FaceAlpha',seaalpha);
 
  yp = y(:);
  yp = [yp(1); yp; yp(end:-1:1)];
  xp = repmat(x(1),size(yp));
  zp = DEM.Z(:,1);
  zp = min(zp(:),0);
  zp = double([sealevel;zp;repmat(sealevel,size(zp))]);
- p(3) = patch(xp,yp,zp,facecolordark,'EdgeColor','none','FaceColor',facecolordark,'FaceAlpha',seaalpha);
+ pa(3) = patch(xp,yp,zp,facecolordark,'EdgeColor','none','FaceColor',facecolordark,'FaceAlpha',seaalpha);
 
  xp = repmat(x(end),size(yp));
  zp = DEM.Z(:,end);
  zp = min(zp(:),0);
  zp = double([sealevel;zp;repmat(sealevel,size(zp))]);
- p(4) = patch(xp,yp,zp,facecolordark,'EdgeColor','none','FaceColor',facecolordark,'FaceAlpha',seaalpha);
+ pa(4) = patch(xp,yp,zp,facecolordark,'EdgeColor','none','FaceColor',facecolordark,'FaceAlpha',seaalpha);
 
  I = DEM.Z>0;
- H = zeros(DEM.size);
- H(I) = nan;
- % Texturemap
- I = ~I;
- TM = repmat(I,1,1,3);
+%  H = zeros(DEM.size);
+%  H(I) = nan;
+%  % Texturemap
+%  I = ~I;
+%  TM = repmat(I,1,1,3);
 
  hold on
  h = surf(x,y,+I,'FaceAlpha',0.5,'FaceColor',facecolor,'EdgeColor','none');
  hold off
 
 
- set(p,'FaceLighting','none')
-
+ set(pa,'FaceLighting','none')
 
+ if ~isempty(baselevel)
  axislim = axis;
  axislim(5) = baselevel;
  axis(axislim)
+ end
 
 
 
@@ -289,7 +280,7 @@ function exaggerate(axes_handle,exagfactor)
 % Author: Wolfgang Schwanghart (w.schwanghart[at]unibas.ch)
 % Date: 6. September, 2010
 
-if nargin == 1;
+if nargin == 1
  exagfactor = 1;
 end
 axis(axes_handle,'image');