hillshade

 PUPROSE: Calculate hillshade for a digital elevation model (DEM)
 -------------------------------------------------------------------
 USAGE: h = hillshade(dem,X,Y,varagin)
 where: dem is the DEM to calculate hillshade for
        X and Y are the DEM coordinate vectors
        varargin are parameters options

 OPTIONS: 
        'azimuth'  is the direction of lighting in deg (default 315)
        'altitude' is the altitude of the lighting source in
                   in degrees above horizontal (default 45)
        'zfactor'  is the DEM altitude scaling z-factor (default 1)
        'plotit'   creates a simple plot of the hillshade

 EXAMPLE:
       h=hillshade(peaks(50),1:50,1:50,'azimuth',45,'altitude',100,'plotit')
       - calculates the hillshade for an example 50x50 peak surface.
       - changes the default settings for azimuth and altitude.
       - creates an output hillshade plot

0001 function h = hillshade(dem,X,Y,varargin)
0002 % PUPROSE: Calculate hillshade for a digital elevation model (DEM)
0003 % -------------------------------------------------------------------
0004 % USAGE: h = hillshade(dem,X,Y,varagin)
0005 % where: dem is the DEM to calculate hillshade for
0006 %        X and Y are the DEM coordinate vectors
0007 %        varargin are parameters options
0008 %
0009 % OPTIONS:
0010 %        'azimuth'  is the direction of lighting in deg (default 315)
0011 %        'altitude' is the altitude of the lighting source in
0012 %                   in degrees above horizontal (default 45)
0013 %        'zfactor'  is the DEM altitude scaling z-factor (default 1)
0014 %        'plotit'   creates a simple plot of the hillshade
0015 %
0016 % EXAMPLE:
0017 %       h=hillshade(peaks(50),1:50,1:50,'azimuth',45,'altitude',100,'plotit')
0018 %       - calculates the hillshade for an example 50x50 peak surface.
0019 %       - changes the default settings for azimuth and altitude.
0020 %       - creates an output hillshade plot
0021 
0022 % See also: GRADIENT, CART2POL
0023 %
0024 % Note: Uses simple unweighted gradient of 4 nearest neighbours for slope
0025 %       calculation (instead of Horn's method) with ESRIs hillshade
0026 %       algorithm.
0027 %
0028 % Felix Hebeler, Dept. of Geography, University Zurich, February 2007.
0029 % modified by Andrew Stevens (astevens@usgs.gov), 5/04/2007
0030 
0031 %% configure inputs
0032 %default parameters
0033 azimuth=315;
0034 altitude=45;
0035 zf=1;
0036 plotit=0;
0037 
0038 %parse inputs
0039 if isempty(varargin)~=1     % check if any arguments are given
0040     [m1,n1]=size(varargin);
0041     opts={'azimuth';'altitude';'zfactor';'plotit'};
0042     for i=1:n1;             % check which parameters are given
0043         indi=strcmpi(varargin{i},opts);
0044         ind=find(indi==1);
0045         if isempty(ind)~=1
0046             switch ind
0047                 case 1
0048                     azimuth=varargin{i+1};
0049                 case 2
0050                     altitude=varargin{i+1};
0051                 case 3
0052                     zf=varargin{i+1};
0053                 case 4
0054                     plotit=1;
0055             end
0056         end
0057     end
0058 end
0059 
0060 %% Initialize paramaters
0061 dx=abs(X(2)-X(1));  % get cell spacing in x and y direction
0062 dy=abs(Y(2)-Y(1));  % from coordinate vectors
0063 
0064 % lighting azimuth
0065 azimuth = 360.0-azimuth+90; %convert to mathematic unit
0066 azimuth(azimuth>=360)=azimuth-360;
0067 azimuth = azimuth * (pi/180); %  convert to radians
0068 
0069 %lighting altitude
0070 altitude = (90-altitude) * (pi/180); % convert to zenith angle in radians
0071 
0072 %% calc slope and aspect (radians)
0073 [fx,fy] = gradient(dem,dx,dy); % uses simple, unweighted gradient of immediate neighbours
0074 [asp,grad]=cart2pol(fy,fx); % convert to carthesian coordinates
0075 %grad = grad/d; % multiply w cellsize
0076 grad=atan(zf*grad); %steepest slope
0077 % convert asp
0078 asp(asp<pi)=asp(asp<pi)+(pi/2);
0079 asp(asp<0)=asp(asp<0)+(2*pi);
0080 
0081 %% hillshade calculation
0082 h = ( (cos(altitude).*cos(grad) ) + ( sin(altitude).*sin(grad).*cos(azimuth-asp)) ); % ESRIs algorithm
0083 h(h<0)=0; % set hillshade values to min of 0.
0084 
0085 h=setborder(h,1,NaN); % set border cells to NaN
0086 
0087 %% plot results if requested
0088 if plotit==1
0089     figure
0090     imagesc(X,Y,h)
0091     axis image
0092     set(gca,'ydir','normal')
0093     colormap(gray)
0094 end
0095 
0096 %% -- Subfunction--------------------------------------------------------------------------
0097 function grid = setborder(grid,bs,bvalue)
0098 grid(1:bs,:)=bvalue; %toprows
0099 grid(size(grid,1)-bs+1:size(grid,1),:)=bvalue; %bottom rows
0100 grid(:,1:bs)=bvalue; %left cols
0101 grid(:,size(grid,2)-bs+1:size(grid,2))=bvalue;