diff --git a/toolbox/@FLOWobj/flipdir.m b/toolbox/@FLOWobj/flipdir.m
index 8227f14..6874f7b 100644
--- a/toolbox/@FLOWobj/flipdir.m
+++ b/toolbox/@FLOWobj/flipdir.m
@@ -74,7 +74,7 @@
 if isempty(options.weights)
     A = flowacc(FD);
     w = A.Z(FD.ix).*FD.fraction;
-elseif ischar(options.weights)
+elseif ischar(options.weights) || isstring(options.weights)
     switch  lower(options.weights)
         case 'rand'
            w = rand(size(FD.fraction));
diff --git a/toolbox/@FLOWobj/reweight.m b/toolbox/@FLOWobj/reweight.m
new file mode 100644
index 0000000..51d36a8
--- /dev/null
+++ b/toolbox/@FLOWobj/reweight.m
@@ -0,0 +1,109 @@
+function FD = reweight(FD,DEM,beta,options)
+
+%REWEIGHT Re-calculate flow weights in a multiple flow direction FLOWobj
+%
+% Syntax
+%
+%     FD = reweight(FD,DEM,beta)
+%     FD = reweight(FD,DEM,beta,pn,pv,...)
+%
+% Description
+%
+%     The function reweight takes a multiple flow direction object FD and
+%     re-calculates the proportions at which flow is distributed from each  
+%     pixel to its downstream neighbors. By default, proportions in FD 
+%     calculated with FLOWobj(DEM,'multi') are proportional to slope, so 
+%     that beta = 1.
+%
+%     Reweight recalculates the slopes based on the DEM and then uses beta
+%     to reweight and normalizes the edges of the flow network. This is the
+%     approach taken by Holmgren (1994).
+%
+%     Note that reweight does not change the topology of the flow network
+%     so that no links are added nor removed. If slopes between two nodes
+%     in the network are zero or negative, their slope is adjusted to have
+%     a very small value (0.1 * the smallest non-negative slope value).
+%
+% Input arguments
+%
+%     FD      FLOWobj (multiple flow directions)
+%     DEM     Digital elevation model (GRIDobj)
+%     beta    exponent >= 0 
+%
+%     Parameter name/value pairs
+%
+%     minslope    scalar (0.1). This indicates the multiplier used to 
+%                 derive the minimum slope value for links that have zero
+%                 or negative slopes.
+%     stable      false or true. Uses a numerically stable approach when
+%                 using large exponents. By default false, but if beta > 5,
+%                 then the stable approach is chosen.
+%     
+% Output arguments
+%
+%     FD      FLOWobj
+%
+% Example
+%
+%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
+%     FD  = FLOWobj(DEM,'multi');
+%     FD1 = reweight(FD,DEM,0);
+%     FD2 = reweight(FD,DEM,10);
+%     A1 = flowacc(FD1);
+%     A2 = flowacc(FD2);  
+%     tiledlayout(1,2,"TileSpacing","compact");
+%     ax(1) = nexttile; imageschs(DEM,log(A1),'colormap',flowcolor);
+%     ax(2) = nexttile; imageschs(DEM,log(A2),'colormap',flowcolor);
+%     linkaxes(ax,'xy');
+%     ext = {[3.7769e+05 3.8353e+05]   [3.7911e+06 3.7948e+06]};
+%     setextent(ext,ax(1))
+%
+% See also: FLOWobj, FLOWobj/multi_normalize
+%  
+% Author: Wolfgang Schwanghart (schwangh[at]uni-potsdam.de)
+% Date: 6. April, 2026
+
+arguments 
+    FD FLOWobj
+    DEM GRIDobj
+    beta (1,1) {mustBeNonnegative}
+    options.minslope (1,1) {mustBeNonnegative} = 0.1
+    options.stable (1,1) = beta > 5
+end
+
+[X,Y] = getcoordinates(DEM,'mat');
+dz = double(DEM.Z(FD.ix)) - double(DEM.Z(FD.ixc));
+dx = sqrt((X(FD.ix) - X(FD.ixc)).^2 + ...
+          (Y(FD.ix) - Y(FD.ixc)).^2);
+FD.fraction = dz./dx;
+minNonZeroFrac = min(FD.fraction(FD.fraction>0));
+FD.fraction = max(FD.fraction,minNonZeroFrac*0.1);
+
+if ~options.stable
+    FD.fraction = FD.fraction.^beta;
+    FD = multi_normalize(FD);
+else
+    % [~,~,locb] = unique(FD.ix);
+    % [W,maxY] = splitapply(@(x) normfun(x,beta),FD.fraction,locb);
+    % FD.fraction = exp(beta.*log(FD.fraction) - maxY(locb))./W(locb);
+    W = accumarray(FD.ix,FD.fraction,[numel(DEM.Z) 1],@(x) normfun(x,beta));
+    maxY = accumarray(FD.ix,FD.fraction,[numel(DEM.Z) 1],@(x) normfun2(x,beta));    
+    FD.fraction = exp(beta.*log(FD.fraction) - maxY(FD.ix))./W(FD.ix);
+
+end
+
+end
+
+function [sumw,maxy] = normfun(x,beta)
+    y = beta * log(x);
+    maxy = max(y);
+    y = y - maxy;
+    w = exp(y);
+    % w = w / sum(w);
+    sumw = sum(w);
+end
+
+function maxy = normfun2(x,beta)
+    y = beta * log(x);
+    maxy = max(y);
+end
\ No newline at end of file
diff --git a/toolbox/@GRIDobj/getcoordinates.m b/toolbox/@GRIDobj/getcoordinates.m
index 853343e..26ee5f4 100644
--- a/toolbox/@GRIDobj/getcoordinates.m
+++ b/toolbox/@GRIDobj/getcoordinates.m
@@ -48,10 +48,6 @@
         [x,y] = meshgrid(x,y);
     case 'GRIDobj'
         [x,y] = meshgrid(x,y);
-        X = GRIDobj(DEM);
-        X.Z = x;
-        Y = GRIDobj(DEM);
-        Y.Z = y;
-        x = X;
-        y = Y;
+        x = GRIDobj(DEM,x);
+        y = GRIDobj(DEM,y);    
 end
\ No newline at end of file
diff --git a/toolbox/@GRIDobj/getoutline.m b/toolbox/@GRIDobj/getoutline.m
index 5a83d07..4cbf29f 100644
--- a/toolbox/@GRIDobj/getoutline.m
+++ b/toolbox/@GRIDobj/getoutline.m
@@ -47,7 +47,7 @@
 %     
 %
 % Author: Wolfgang Schwanghart (schwangh[at]uni-potsdam.de)
-% Date: 30. May, 2024
+% Date: 20. April, 2026
 
 
 arguments (Input)
@@ -61,13 +61,14 @@
 nnan = options.shownans && any(I);
 
 if ~nnan
-    csh   = DEM.cellsize/2;
-    [x,y] = getcoordinates(DEM);
-    maxx = max(x)+csh/2;
-    minx = min(x)-csh/2;
     
-    maxy = max(y)+csh/2;
-    miny = min(y)-csh/2;
+    [x,y] = getcoordinates(DEM);
+
+    csh   = DEM.cellsize/2;
+    maxx = max(x)+csh;
+    minx = min(x)-csh;    
+    maxy = max(y)+csh;
+    miny = min(y)-csh;
     
     x = [minx minx maxx maxx minx];
     y = [miny maxy maxy miny miny];
@@ -79,15 +80,19 @@
         GT.Properties.VariableNames = {'Lat','Lon'};
         GT = table2geotable(GT,"CoordinateReferenceSystem",DEM.georef.GeographicCRS,...
             "GeometryType","polygon");
+        GT = removevars(GT,["Lat" "Lon"]);
+        
     elseif isProjected(DEM)
         GT = table(x,y);
         GT.Properties.VariableNames = {'X','Y'};
         GT = table2geotable(GT,"CoordinateReferenceSystem",DEM.georef.ProjectedCRS,...
             "GeometryType","polygon");
+        GT = removevars(GT,["X" "Y"]);
     else
         GT = table(x,y);
         GT.Properties.VariableNames = {'X','Y'};
         GT = table2geotable(GT,"GeometryType","polygon");
+        GT = removevars(GT,["X" "Y"]);
 
     end
     else
diff --git a/toolbox/@GRIDobj/inpaintnans.m b/toolbox/@GRIDobj/inpaintnans.m
index 1f209d1..349797c 100644
--- a/toolbox/@GRIDobj/inpaintnans.m
+++ b/toolbox/@GRIDobj/inpaintnans.m
@@ -331,32 +331,33 @@
 %     DEMi     GRIDobj with
 %
 
-
-I = isnan(DEM.Z);
-sq2 = sqrt(2);
-w = 1./[sq2 1 sq2; ...
-     1   0   1; ...
-     sq2 1 sq2];
-w = w(:);
-
-Z = nlfilter(DEM.Z,[3 3],@fun);
-DEM.Z = Z;
-
-function b = fun(a)
-
-a = a(:);    
-if ~isnan(a(5))
-    b = a(5);
-    return
-end
-
-I = isnan(a);
-if all(I)
-    b = nan;
-    return
-end
-
-I = ~I;
-b = sum(a(I).*(w(I)./sum(w(I))));
-end
+DEM.Z = fillmissing2(DEM.Z,"movmean",3);
+
+% I = isnan(DEM.Z);
+% sq2 = sqrt(2);
+% w = 1./[sq2 1 sq2; ...
+%      1   0   1; ...
+%      sq2 1 sq2];
+% w = w(:);
+% 
+% Z = nlfilter(DEM.Z,[3 3],@fun);
+% DEM.Z = Z;
+% 
+% function b = fun(a)
+% 
+% a = a(:);    
+% if ~isnan(a(5))
+%     b = a(5);
+%     return
+% end
+% 
+% I = isnan(a);
+% if all(I)
+%     b = nan;
+%     return
+% end
+% 
+% I = ~I;
+% b = sum(a(I).*(w(I)./sum(w(I))));
+% end
 end
diff --git a/toolbox/@GRIDobj/pad.m b/toolbox/@GRIDobj/pad.m
index ce34a07..6057c90 100644
--- a/toolbox/@GRIDobj/pad.m
+++ b/toolbox/@GRIDobj/pad.m
@@ -20,8 +20,10 @@
 % Input arguments
 %
 %     DEM     GRIDobj
-%     px      amount of padding on each of the four edges of the grid. 
-%             For px=1 the resulting size of DEMp is DEM.size+2. 
+%     px      scalar or 1x4 vector indicating the amount of padding on each 
+%             of the four edges of the grid. For px=1 the resulting size of
+%             DEMp is DEM.size+2. If px is a vector, then elements indicate
+%             the amount of padding in [north south west east] direction.             
 %     val     pad value (default=0). This value does not have any effects 
 %             if px<0.
 %
@@ -46,18 +48,22 @@
 
 arguments
     DEM   GRIDobj
-    px    {mustBeNumeric} = 1
+    px    (1,4) {mustBeNumeric} = 1
     val   {mustBeScalarOrEmpty} = 0
 end
 
-if px == 0
+% If all px are zero, we can return the original GRIDobj
+if all(px == 0)
     return
 end
 
-px  = sign(px)*ceil(abs(px));
+% Make sure that we have integer values
+px  = sign(px).*ceil(abs(px));
 
 % new size of output grid
-newsize = DEM.size + px*2;
+% px(1): north px(2):south px(3): west px(4): east
+newsize = [DEM.size(1) + px(1) + px(2) ... 
+           DEM.size(2) + px(3) + px(4)];
 
 % check if size is negative and issue an error
 if any(newsize < 0)
@@ -75,16 +81,22 @@
 end
 
 % transfer values to new array
-if px>0
-    Znew(px+1:end-px,px+1:end-px) = DEM.Z;
-elseif px<0
+if all(px>0)
+    Znew(px(1)+1:end-px(2),px(3)+1:end-px(4)) = DEM.Z;
+elseif all(px<0)
     abspx = abs(px);
-    Znew = DEM.Z(abspx+1:end-abspx,abspx+1:end-abspx);
+    Znew = DEM.Z(abspx(1)+1:end-abspx(2),abspx(3)+1:end-abspx(4));
+else
+    pxpos = max(px,0);
+    pxneg = -min(px,0);
+
+    Znew(pxpos(1)+1:end-pxpos(2),pxpos(3)+1:end-pxpos(4)) = ...
+        DEM.Z(pxneg(1)+1:end-pxneg(2),pxneg(3)+1:end-pxneg(4));
 end
 
 % adjust referencing
 DEM.Z           = Znew;
-DEM.wf(:,3)     = DEM.wf(:,3) -px.*[DEM.wf(1,1) DEM.wf(2,2)]';
+DEM.wf(:,3)     = DEM.wf(:,3) -[px(3) px(1)]'.*[DEM.wf(1,1) DEM.wf(2,2)]';
 DEM.size        = size(DEM.Z);
 
 if ~isempty(DEM.georef)
@@ -106,7 +118,7 @@
             case 'map.rasterref.MapCellsReference'
                 Rnew = maprasterref(DEM.wf,DEM.size,R.RasterInterpretation);
             case 'map.rasterref.GeographicCellsReference'
-            Rnew = georasterref(DEM.wf,DEM.size,R.RasterInterpretation);
+                Rnew = georasterref(DEM.wf,DEM.size,R.RasterInterpretation);
         end
     end
 
diff --git a/toolbox/@PPS/PPS.m b/toolbox/@PPS/PPS.m
index 11cef0f..1654a99 100644
--- a/toolbox/@PPS/PPS.m
+++ b/toolbox/@PPS/PPS.m
@@ -259,7 +259,7 @@
             elseif isa(PS,'geoshape')
                 lat = PS.Latitude;
                 lon = PS.Longitude;
-                [X1,Y1] = mfwdtran(P.S.georef.mstruct,lat,lon);
+                [X1,Y1] = projfwd(P.S.georef.ProjectedCRS,lat,lon);
             end
             [X,Y] = STREAMobj2XY(P.S);
             [Xi,Yi] = polyxpoly(X1,Y1,X,Y);
diff --git a/toolbox/@STREAMobj/crs.m b/toolbox/@STREAMobj/crs.m
index 64da846..f7f3b24 100644
--- a/toolbox/@STREAMobj/crs.m
+++ b/toolbox/@STREAMobj/crs.m
@@ -187,6 +187,9 @@
     % a special case. The stream network needs at least three nodes in a row. If
     % there are less, the second derivative matrix is empty, and crs uses
     % quantile carving instead
+    if isnan(options.mingradient)
+        options.mingradient = 0;
+    end
     zs = quantcarve(S,z,options.tau,...
                         'split',options.split,...
                         'mingradient',options.mingradient,...