several updates

@GRIDobj/surf.m

@@ -80,7 +80,9 @@
 if block
  minz = min(DEM);
  maxz = max(DEM);
+ if isempty(baselevel)
  baselevel = minz-(maxz-minz)*0.2;
+ end
 end 
 sea = p.Results.sea;
 sealevel = p.Results.sealevel;

@PPS/PPS.m

@@ -8,6 +8,7 @@
 % P = PPS(S,'PP',xy,...)
 % P = PPS(S,'PP',MS,...)
 % P = PPS(S,'PP',nal,...)
+% P = PPS(S,'PP',geotable,...)
 % P = PPS(S,'runif',n,...)
 % P = PPS(S,'rpois',lambda,...)
 % P = PPS(S,'intersect',PS,...)
@@ -107,14 +108,14 @@
 % 
 % Reference: Schwanghart, W., Molkenthin, C., & Scherler, D. (2020). A 
 % systematic approach and software for the analysis of point patterns on 
-% river networks. Earth Surface Processes and Landforms, accepted. 
+% river networks. Earth Surface Processes and Landforms, 46 (9), 1847-1862. 
 % [DOI: 10.1002/esp.5127]
 %
 % See also: GRIDobj, FLOWobj, STREAMobj, STREAMobj/snap2stream, SWATHobj
 % DIVIDEobj
 %
 % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
-% Date: 22. January, 2020
+% Date: 21. December, 2021
 
 properties
 %S instance of STREAMobj
@@ -201,6 +202,17 @@
  else
  P.z = results.z;
  end
+
+ %% Enable reading of geotable
+ if ~verLessThan('map','5.2')
+ if isgeotable(results.PP)
+ points = results.PP.Shape;
+ x = [points.X];
+ y = [points.Y];
+
+ results.PP = [x(:) y(:)];
+ end
+ end
 
  %% Read or create points
  if ~isempty(results.PP) 

@PPS/aggregate.m

@@ -1,4 +1,4 @@
-function [P,locb,n] = aggregate(P,c,varargin)
+function [P,locb] = aggregate(P,c,varargin)
 
 %AGGREGATE Aggregate points in PPS to new point pattern
 %

@PPS/fitloglinear.m

@@ -23,7 +23,7 @@
 % 'stepwise' {false} or true. If true, fitloglinear uses stepwiseglm
 % to fit the model.
 % 'modelspec' see fitglm. For example, for fitting a forth-order
-% polynomial: 'poly4'
+% polynomial of one covariate: 'poly4'
 % 
 % In addition, fitloglinear accepts parameter name/value pairs of
 % fitglm or stepwiseglm.

@PPS/roc.m

@@ -78,10 +78,12 @@
 pp = points(P,'nal');
 
 if p.Results.perfcurve
+ TVals = quantile(c,linspace(0,1,p.Results.nTvals));
+ % TVals = linspace(min(c),max(c),p.Results.nTvals);
  [X,Y,~,AUC] = perfcurve(pp,c,true,...
  'NBoot',p.Results.NBoot,...
  'BootType',p.Results.BootType,...
- 'TVals',linspace(min(c),max(c),p.Results.nTvals));
+ 'TVals',TVals);
 else
 
  cp = sortrows([c pp],[1 2]);

@STREAMobj/ksn.m

@@ -83,7 +83,7 @@
 end
 
 % minima imposition to avoid negative gradients
-z = imposemin(S,z);
+z = imposemin(S,z,0.00001);
 % calculate gradient
 g = gradient(S,z);
 % upslope area
@@ -94,7 +94,10 @@
 k = g./(a.^(-p.Results.theta));
 
 if p.Results.smooth ~= 0
- k = smooth(S,k,'K',p.Results.smooth);
+ k(k==0) = 0.0001;
+ ks = smooth(S,log(k),'K',p.Results.smooth);
+% sig = var(ks-log(k));
+ k = exp(ks)*exp(var(-ks/2));
 end
 
 

@STREAMobj/mchi.m

@@ -11,8 +11,8 @@
 % mchi calculates the gradient of elevation in a DEM where the
 % horizontal coordinate is chi as obtained from the function
 % chitransform. mchi provides a metric that can be used to compare
-% channel gradients with different drainage areas and is very similar
-% to ksn.
+% channel gradients with different drainage areas and is the same as
+% ksn.
 %
 % Input arguments
 %

IOtools/readopentopo.m

@@ -15,6 +15,10 @@
 % projected to a projected coordinate system (use reproject2utm) before
 % analysis in TopoToolbox.
 %
+% NOTE: Starting on January 1st, 2022, an API authorization key will be 
+% required for this API. Users can request an API key via myOpenTopo in 
+% the OpenTopography portal (https://opentopography.org/developers).
+%
 % Input arguments
 %
 % Parameter name values
@@ -36,17 +40,17 @@
 % 'south' southern boundary
 % 'west' western boundary
 % 'east' eastern boundary
-% 'demtype' The global raster dataset 
-% {'SRTMGL3'}: SRTM GL3 (90m) (default) 
+% 'demtype' The global raster dataset *
+% {'SRTMGL3'}: SRTM GL3 (90m) (default)
 % 'SRTMGL1': SRTM GL1 (30m) 
 % 'SRTMGL1_E': SRTM GL1 (Ellipsoidal) 
 % 'AW3D30': ALOS World 3D 30m 
 % 'AW3D30_E': ALOS World 3D (Ellipsoidal)
 % 'SRTM15Plus': Global Bathymetry SRTM15+ V2.1 (only 
 % mediterranean area so far)
-% 'NASADEM': NASADEM Global DEM *
-% 'COP30': Copernicus Global DSM 30m *
-% 'COP90': Copernicus Global DSM 90m *
+% 'NASADEM': NASADEM Global DEM 
+% 'COP30': Copernicus Global DSM 30m 
+% 'COP90': Copernicus Global DSM 90m 
 % 
 % * requires API Key (see option 'apikey').
 %
@@ -110,25 +114,23 @@
 f = fullfile(p.Results.filename);
 
 % check api
-if any(strcmp(demtype,{'NASADEM','COP30','COP90'}))
- if isempty(p.Results.apikey)
- % check whether file opentopography.apikey is available
- if exist('opentopography.apikey','file')
- fid = fopen('opentopography.apikey');
- apikey = textscan(fid,'%c');
- apikey = apikey{1}';
- % Remove trailing blanks, if there are any
- apikey = deblank(apikey);
- else
- error('The DEM types NASADEM, COP30 and COP90 require an API Key')
- end
+
+if isempty(p.Results.apikey)
+ % check whether file opentopography.apikey is available
+ if exist('opentopography.apikey','file')
+ fid = fopen('opentopography.apikey');
+ apikey = textscan(fid,'%c');
+ apikey = apikey{1}';
+ % Remove trailing blanks, if there are any
+ apikey = deblank(apikey);
  else
- apikey = p.Results.apikey;
+ error('The DEM types NASADEM, COP30 and COP90 require an API Key')
  end
 else
- apikey = [];
+ apikey = p.Results.apikey;
 end
 
+
 % save to drive
 options = weboptions('Timeout',100000);