directsun

 PURPOSE: Calculate binary grid for a digital elevation model (DEM) static
          whether a cell is shaded or in direct sunlight.
 -------------------------------------------------------------------
 USAGE: h = directsun(dem,X,Y,azimuth,altitude,sf)
 where: dem is the DEM to calculate radiation for
        X and Y are the coordinate vectors
 Options: provide via varargin ('option',value,...)
        'azimuth' is the direction of lighting source in degrees (default 315deg)
        'altitude' is the altitude of the lighting source in
        in degrees above horizontal (default 45deg)
        'sf' is the shading factor between 0..1 deciding what is
        considered shading
        'zf' is the z-factor    


 Felix Hebeler, Geography Dept., University Zurich, May 2007.

0001 function h = directsun(dem,X,Y,varargin)
0002 % PURPOSE: Calculate binary grid for a digital elevation model (DEM) static
0003 %          whether a cell is shaded or in direct sunlight.
0004 % -------------------------------------------------------------------
0005 % USAGE: h = directsun(dem,X,Y,azimuth,altitude,sf)
0006 % where: dem is the DEM to calculate radiation for
0007 %        X and Y are the coordinate vectors
0008 % Options: provide via varargin ('option',value,...)
0009 %        'azimuth' is the direction of lighting source in degrees (default 315deg)
0010 %        'altitude' is the altitude of the lighting source in
0011 %        in degrees above horizontal (default 45deg)
0012 %        'sf' is the shading factor between 0..1 deciding what is
0013 %        considered shading
0014 %        'zf' is the z-factor
0015 %
0016 %
0017 % Felix Hebeler, Geography Dept., University Zurich, May 2007.
0018 
0019 sf=0.3; % shadingfactor: set to value between 0..1 to determine what is considered
0020 azimuth=315;
0021 altitude=45;
0022 zf=1;
0023 
0024 %% inputs
0025 if isempty(varargin)~=1     % check if any arguments are given
0026     [m1,n1]=size(varargin);
0027     opts={'azimuth';'altitude';'sf','zf'};
0028     for i=1:n1;             % check which parameters are given
0029         indi=strcmpi(varargin{i},opts);
0030         ind=find(indi==1);
0031         if isempty(ind)~=1
0032             switch ind
0033                 case 1
0034                     azimuth=varargin{i+1};
0035                 case 2
0036                     altitude=varargin{i+1};
0037                 case 3
0038                     sf=varargin{i+1};
0039                 case 4
0040                     zf=varargin{i+1};
0041 
0042             end
0043         end
0044     end
0045 end
0046 
0047 %% Initialize paramaters
0048 %h = nan(size(dem,1),size(dem,2));
0049 
0050 dx=abs(X(2)-X(1));
0051 dy=abs(Y(2)-Y(1));
0052 
0053 % lighting azimuth
0054 azimuth = 360.0-azimuth+90; %convert to mathematic unit
0055 azimuth(azimuth>=360)=azimuth-360;
0056 azimuth = azimuth * (pi/180); %  convert to radians
0057 
0058 %lighting altitude
0059 altitude = (90-altitude) * (pi/180); % convert to zenith angle and radians
0060 
0061 %% calc slope and aspect (radians)
0062 [fx,fy] = gradient(dem,dx,dy); % uses simple, unweighted gradient of immediate neighbours
0063 [asp,grad]=cart2pol(fy,fx); % convert to carthesian coordinates
0064 %grad = grad/d; % multiply w cellsize
0065 grad=atan(zf*grad); %steepest slope
0066 % convert asp
0067 asp(asp<pi)=asp(asp<pi)+(pi/2);
0068 asp(asp<0)=asp(asp<0)+(2*pi);
0069 
0070 %% hillshade calculation
0071 h = ( (cos(altitude).*cos(grad) ) + ( sin(altitude).*sin(grad).*cos(azimuth-asp)) ); % ESRIs algorithm
0072 h(h<0)=0;
0073 h(h<sf)=0;
0074 h(h>=sf)=1;
0075 h=setborder(h,1,NaN);
0076 
0077 %--------------------------------------------------------------------------
0078 function grid = setborder(grid,bs,bvalue)
0079 grid(1:bs,:)=bvalue; %toprows
0080 grid(size(grid,1)-bs+1:size(grid,1),:)=bvalue; %bottom rows
0081 grid(:,1:bs)=bvalue; %left cols
0082 grid(:,size(grid,2)-bs+1:size(grid,2))=bvalue;