LS code in matlab from rafa

LS_from_DEM.m

clear all
close all

addpath(genpath(pwd))

% channel init threshold [km2]
Ad_ini=1;

% length threshold [m]

Lcap=333;

% load input data that are already gridded
DEM=GRIDobj('DEM_hydrosheds_CR_clip_watershed_197547.tif');

% get flow directions, accumulations, distances
FDir=FLOWobj(DEM);
FDirInf=FLOWobj(DEM,'Dinf');

Ad=flowacc(FDir);
%GRIDobj2geotiff(Ad,[pwd '/Outputs_ML/' 'AD_cells.tif'])
Ad=Ad.*DEM.cellsize.^2;
%GRIDobj2geotiff(Ad,[pwd '/Outputs_ML/' 'AD_m2.tif'])

% START THE USLE PROCEDURE
AdDInf=flowacc(FDirInf).*DEM.cellsize.^2;
fDist = flowdistance(FDir,'downstream');

% get gradients
slope_perc=gradient8(DEM,'perc'); % note that many of the slope filtering for SDR requires degree values
theta=gradient8(DEM,'deg'); % most other calculations need either degree or radian values. Here I use degree. Note that we then need to use "sind()" instead of "sin()".

% calculate m
m=GRIDobj(DEM); % initialize storage

% SDR approach
    ii=slope_perc.Z<=1;
        m.Z(ii)=0.2;

    ii=slope_perc.Z > 1 & slope_perc.Z<=3.5;
        m.Z(ii)=0.3;

    ii=slope_perc.Z > 3.5 & slope_perc.Z<=5;
        m.Z(ii)=0.4;

    ii=slope_perc.Z >  5 & slope_perc.Z<=9;
        m.Z(ii)=0.5;

    ii=slope_perc.Z>9;
        beta=sind(theta.Z)./0.0986./(3.*sind(theta.Z.^0.8)+0.56);

% Borelli approach. Results in really large values
%         beta=sind(theta.Z)./0.0986./(3.*sind(theta.Z.^0.8)+0.56);
%         m.Z=beta./(1+beta);

% calculate L factor
    % inputs
    D=DEM.cellsize;
    a=aspect(DEM,'deg');
    a.Z=double(a.Z);
    x=sind(a.Z)+cosd(a.Z);

    %Borrelli approach
        Ain=Ad;
        Ain.Z(fDist.Z>Lcap)=Lcap.^2; % note conversion between areas and slopes using the square root.
        L=((Ain+D.^2)^(m+1)-Ain.^(m+1))./(D.^(m+2).*x.^m.*22.13.^m);

        % Invest
        Ain_inv=Ad;
        L_inv=((Ain_inv+D.^2)^(m+1)-Ain_inv.^(m+1))./(D.^(m+2).*x.^m.*22.13.^m);
        L_inv.Z(L_inv.Z>Lcap)=Lcap;

    % calculate S factor
        S=GRIDobj(DEM); % initializs S

        ii=slope_perc.Z<9;
        S.Z(ii)=10.8*sind(theta.Z(ii))+0.03;

        ii=slope_perc.Z>=9;
        S.Z(ii)=16.8*sind(theta.Z(ii))-0.5;

%     figure; histogram(slope_perc.Z(slope_perc.Z>0))
%     figure; histogram(S.Z(S.Z>0))

% calculate LS factor
LS=L*S;
LS.Z(LS.Z==0)=nan

LS_inv=L_inv*S;
LS_inv.Z(LS_inv.Z==0)=nan

figure
subplot(1,2,1)
imageschs(LS)
subplot(1,2,2)
imageschs(LS_inv)

GRIDobj2geotiff(S, 'S_from_ML.tif')
GRIDobj2geotiff(LS, 'LS_from_ML.tif')
GRIDobj2geotiff(LS, 'LS_from_ML_InVEST_pre3_7_style.tif')