DEM_Detrend.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------DEM Detrender 1.0 -----------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Produced by Joseph Wheaton                 %
%                           July 2007                        %
%                                                                %
%               Last Updated: 13 July, 2007                  %
%                                                                %
%                                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This program detrends the elvation values in a reach DEM based on
% calculated valley slope. Assumes grid is rotated such that valley
% is sloped from high points in south to low points in north.
%


% Edit History

%------------------------------------------------------------------- 
	clc
    clear all
    echo off
	tic                                             % Start Timer
%------Read in Data -------------------------------------------------
       h = msgbox('This program detrends the elvation values in a reach DEM based on specified valley slope. NOTE: the valley slope parameter must be specified in the code.','DEM Detrend','help');
       uiwait(h);
 
    % Read first grid data
    
    [filename, pathname]=uigetfile('*.asc','Select the DEM to detrend (ARC ascii format)');    
    grid1_filename=[pathname filename];
    fid=fopen(grid1_filename,'r');      % opens file to fid in read only mode
    dum1=fscanf(fid,'%s',1);            % assigns first part of header (ncols) to dummy variable dum1 (%s is string notation)
    nx=fscanf(fid,'%u',1);              % stores actual number of collumns to nx (%u is decimal notation)
    dum2=fscanf(fid,'%s',1);            % assigns second line header info (nrows) to dummy variable dum2
    ny=fscanf(fid,'%u',1);              % stores actual number of rows to ny
    dum3=fscanf(fid,'%s',1);            % assigns third line header info (xllcorner) to dummy variable dum3
    xll=fscanf(fid,'%f',1);             % stores actual lower left x coordinate corner to xll (%f is fixed point notation)
    dum4=fscanf(fid,'%s',1);            % assigns fourth line header info (yllcorner) to dummy variable dum4
    yll=fscanf(fid,'%f',1);             % stores actual lower left y coordinate corner to yll
    dum5=fscanf(fid,'%s',1);            % assigns fifth line header info (cellsize) to dummy variable dum5
    lx=fscanf(fid,'%f',1);              % stores cell size in lx
    dum6=fscanf(fid,'%s',1);            % assigns sixth line header info (nodata) to dummy variable dum6
    nodata=fscanf(fid,'%f',1);          % stores the no data ARC tag to nodata
    grid1=fscanf(fid,'%f',[nx,ny]);      % stores all of the cell data in a double array dimensioned according to nx and ny
    fclose(fid);                        % closes the file
    fprintf('Done reading ARC data: %s \n',filename);
    fprintf('\n')
         

    %------Set No-Data to NaN-------------------------------------
    
    nd_cells=find(grid1 == nodata); % Find nodata cells
    grid1(nd_cells)=nan;
    
    % Make new DEM
    dtDEM = zeros(nx,ny);
    dtDEM(nd_cells) = nan;
    
    
    % Find relief for later error checking
    relief = max(max(grid1))- min(min(grid1));
    
%------Parametrs -------------------------------------------------

    
% SPECIFY VALLEY SLOPE MANUALLY!
slope = 0.009;


%------Build Detrended DEM, collumn by collumn--------------------------
    i = 0;
    j = 0;
    rowPosCen = ny/2;
    rowPosBot = 0;
    rowPosTop = 0;
    
    for j=1:nx; % Collumn loop
        % GET START AND STOP POSITIONS
        % First stop from top
        for i=1:ny; % Row loop going down
            if(grid1(j,i) == NaN);
                % Do nothing... advance to next row
            else
                rowPosTop = i;
                break;
            end
        end
        
        % Next stop from bottom
        i=ny;
        for i = ny; -1; 0; % Row loop going up
            if(grid1(j,i) == NaN);
                % Do nothing... advance to next row
            else
                rowPosBot = i;
                break;
            end   
            
        end
        
        if(rowPosBot == 0 | rowPosTop ==0)
            % SKIP ROW
            break;
        else

            % Work from top (with elevation lower)
            i=0;
            for i=rowPosTop:rowPosBot; % Row loop going down
                if(grid1(j,i) == NaN);
                    dtDEM(j,i) == NaN;
                else
                    % Calculate distance of current cell from centre
                    dist = (rowPosCen - i)*lx;
                    % Calculate amount to add or subtract
                    detrend = dist * slope;
                    if(detrend > 0 & detrend > relief | detrend < 0 & abs(detrend) > relief);
                        fprintf('WARNING! At x = %d, y = %d, the elevation adjustment (%5.2f) was greater than the relief (%5.2f)! This cell was not adjusted!',j,i,detrend,relief);
                        detrend = 0; % Reset value and warn user
                    end
                    % Adjust Elevation Values
                    dtDEM(j,i) = grid1(j,i) + detrend;
                end
            end
            
        end
    end


%------Save data to new files ---------------------------------------------------    
    
   
    [filename,pathname]=uiputfile('*.asc','Specify the output file for the detrended DEM');    % Select final file name
    Output_file_name=[pathname filename];
 
    % Switch NaNs back to NoData
    dtDEM(nd_cells) = nodata;
    
    %Write grid1 data to a temporary file in ARC format
    fid1= fopen('temp.txt', 'w');
        i=0;
        j=0;
        for i=1:ny;              % Begin main gridcell loop (rows)
            for j=1:nx;          % Begin gridcell loop (collumns)
                if ((dtDEM(j,i) == nodata) | (dtDEM(j,i) == 0));
                    fprintf(fid1, '%d ',dtDEM(j,i));  
                else    
                    fprintf(fid1, '%6.4f ',dtDEM(j,i));   
                end
            end                 % End gridcell loop (collumns)
            fprintf(fid1, '\n');
        end                     % End main gridcell loop (rows)
    fclose(fid1);
 
    fid1 = fopen('temp.txt', 'r');
    fid2 = fopen(Output_file_name, 'w');
    
    %% Write header in final file
    fprintf(fid2,'ncols \t %d \n', nx);
    fprintf(fid2,'nrows \t %d \n', ny);
    fprintf(fid2,'xllcorner \t %10.4f \n', xll);
    fprintf(fid2,'yllcorner \t %10.4f \n', yll);
    fprintf(fid2,'cellsize \t %d \n', lx);
    fprintf(fid2,'NODATA_value \t %d \n', nodata);
    
    while ~feof(fid1),
        st = fgetl(fid1);
        fprintf(fid2,'%s\n',st);
    end
    fclose(fid1); fclose(fid2);
    delete ('temp.txt');
    

    
     
   

        

%------Clean up and close ---------------------------------------------------
       
        elapsed=toc;
        fprintf('Done writing file (%s). Total time: %5.2f [seconds]\n',filename,elapsed);