xyaspect

 PURPOSE:  computes the polar coordinate x and y components of aspect for a grid
 -------------------------------------------------------------------
 USAGE: [x,y]=xyaspect(Z)
 where: [Z] is a matrix (input DEM) to calculate aspect for.
 RETURNS:
        [X,Y] = are two matrices same size as input matrix Z with the x
        and y coordinates of aspect.
 -------------------------------------------------------------------
 NOTES: X=1 equals East, -1 West
        Y=1 equals North, -1 South
            so X=1/Y=1 is NE, etc

 Felix Hebeler, Geography Dept., University Zurich, July; 2006.

0001 function [X,Y]=xyaspect(Z)
0002 % PURPOSE:  computes the polar coordinate x and y components of aspect for a grid
0003 % -------------------------------------------------------------------
0004 % USAGE: [x,y]=xyaspect(Z)
0005 % where: [Z] is a matrix (input DEM) to calculate aspect for.
0006 % RETURNS:
0007 %        [X,Y] = are two matrices same size as input matrix Z with the x
0008 %        and y coordinates of aspect.
0009 % -------------------------------------------------------------------
0010 % NOTES: X=1 equals East, -1 West
0011 %        Y=1 equals North, -1 South
0012 %            so X=1/Y=1 is NE, etc
0013 %
0014 % Felix Hebeler, Geography Dept., University Zurich, July; 2006.
0015 
0016 X=nlfilter(Z,[3 3],@fxaspect);
0017 Y=nlfilter(Z,[3 3],@fyaspect);
0018 %---------------
0019 function [x] =fxaspect(z)
0020     dsdx = sum(z(1:3))-sum(z(7:9));
0021     dsdy = z(1)+z(4)+z(7)-z(3)-z(6)-z(9);
0022     dsdx2=dsdx*dsdx;
0023     dsdy2=dsdy*dsdy;
0024     if (dsdx2+dsdy2)==0;
0025         x=NaN;
0026         return
0027     end
0028     x=dsdx/sqrt(dsdx2+dsdy2);
0029 
0030     %---------------
0031 function [y] =fyaspect(z)
0032     dsdx = sum(z(1:3))-sum(z(7:9));
0033     dsdy = z(1)+z(4)+z(7)-z(3)-z(6)-z(9);
0034     dsdx2=dsdx*dsdx;
0035     dsdy2=dsdy*dsdy;
0036     if (dsdx2+dsdy2)==0;
0037         y=NaN;
0038         return
0039     end
0040     y=dsdy/sqrt(dsdx2+dsdy2);