DeriveAspect Subroutine

public subroutine DeriveAspect(dem, aspect)

Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbors. Aspect can be thought of as the slope direction. The values of the output raster will be the compass direction of the aspect. A moving window visits each cell in the input raster and for each cell in the center of the window, an aspect value is calculated using an algorithm that incorporates the values of the cell's eight neighbors. If any neighborhood cells are NoData, they are first assigned the value of the center cell, then the aspect is computed. The cells are identified as letters 'a' to 'i', with 'e' representing the cell for which the aspect is being calculated. The rate of change in the x direction for cell 'e' is calculated with the following algorithm:

    [dz/dx] = ((c + 2f + i) - (a + 2d + g)) / 8

The rate of change in the y direction for cell 'e' is calculated with the following algorithm:

    [dz/dy] = ((g + 2h + i) - (a + 2b + c))  / 8

Taking the rate of change in both the x and y direction for cell 'e', aspect is calculated using:

    aspect = 57.29578 * atan2 ([dz/dy], -[dz/dx])

The aspect value is then converted to compass direction values (0–2pi rad), according to the following rule:

    if aspect < 0

    cell = pi/2 - aspect

    else if aspect > pi/2

    cell = 2pi - aspect + pi/2

    else

    cell = pi/2 - aspect

Flat areas in the input raster—areas where the slope is zero are assigned an aspect of -1.

Reference:

Burrough, P. A. and McDonell, R.A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), p. 190

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem
type(grid_real), intent(out) :: aspect

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: a

elevation of N-W cell
real(kind=float), public :: aspectIJ
real(kind=float), public :: b

elevation of N cell
real(kind=float), public :: c

elevation of N-E cell
real(kind=float), public :: d

elevation of W cell
real(kind=float), public :: dZdX

rate of change in the x direction
real(kind=float), public :: dZdY

rate of change in the y direction
real(kind=float), public :: f

elevation of E cell
real(kind=float), public :: g

elevation of S-W cell
real(kind=float), public :: h

elevation of S cell
real(kind=float), public :: i

elevation of S-E cell
integer, public :: ii
integer, public :: jj

Source Code

SUBROUTINE DeriveAspect &
!
(dem, aspect)


IMPLICIT NONE

!arguments with intent in
TYPE(grid_real),    INTENT(in):: dem

!arguments with intent out
TYPE (grid_real), INTENT (out) :: aspect !(rad)

!local variables
INTEGER :: ii,jj
REAL (KIND = float) :: dZdX !!rate of change in the x direction
REAL (KIND = float) :: dZdY !!rate of change in the y direction
REAL (KIND = float) :: aspectIJ
REAL (KIND = float) :: a !!elevation of N-W cell
REAL (KIND = float) :: b !!elevation of N cell
REAL (KIND = float) :: c !!elevation of N-E cell
REAL (KIND = float) :: d !!elevation of W cell
REAL (KIND = float) :: f !!elevation of E cell
REAL (KIND = float) :: g !!elevation of S-W cell
REAL (KIND = float) :: h !!elevation of S cell
REAL (KIND = float) :: i !!elevation of S-E cell

!------------------------------end of declaration -----------------------------
!allocate aspect grid
CALL NewGrid (aspect,dem)

!loop through dem grid
idim: DO ii = 1, dem % idim
  jdim: DO jj = 1, dem % jdim
    IF (dem % mat(ii,jj) /= dem % nodata) THEN
        
        !north-west cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj-1,dem) ) THEN
          IF ( dem % mat (ii-1,jj-1) /= dem % nodata) THEN
            a = dem % mat (ii-1,jj-1)
          ELSE
            a = dem % mat (ii,jj)
          END IF
        ELSE  
          a = dem % mat (ii,jj)  
        END IF
        
        !north cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj,dem) ) THEN
          IF ( dem % mat (ii-1,jj) /= dem % nodata) THEN
            b = dem % mat (ii-1,jj)
          ELSE
            b = dem % mat (ii,jj)
          END IF
        ELSE  
          b = dem % mat (ii,jj)  
        END IF
        
        !north-east cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj+1,dem) ) THEN
          IF ( dem % mat (ii-1,jj+1) /= dem % nodata) THEN
            c = dem % mat (ii-1,jj+1)
          ELSE
            c = dem % mat (ii,jj)
          END IF
        ELSE  
          c = dem % mat (ii,jj)  
        END IF
        
        !west cell
        IF ( .NOT. IsOutOfGrid (ii,jj-1,dem) ) THEN
          IF ( dem % mat (ii,jj-1) /= dem % nodata) THEN
            d = dem % mat (ii,jj-1)
          ELSE
            d = dem % mat (ii,jj)
          END IF
        ELSE  
          d = dem % mat (ii,jj)  
        END IF
        
        !east cell
        IF ( .NOT. IsOutOfGrid (ii,jj+1,dem) ) THEN
          IF ( dem % mat (ii,jj+1) /= dem % nodata) THEN
            f = dem % mat (ii,jj+1)
          ELSE
            f = dem % mat (ii,jj)
          END IF
        ELSE  
          f = dem % mat (ii,jj)  
        END IF
        
        !south-west cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj-1,dem) ) THEN
          IF ( dem % mat (ii+1,jj-1) /= dem % nodata) THEN
            g = dem % mat (ii+1,jj-1)
          ELSE
            g = dem % mat (ii,jj)
          END IF
        ELSE  
          g = dem % mat (ii,jj)  
        END IF
        
        !south cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj,dem) ) THEN
          IF ( dem % mat (ii+1,jj) /= dem % nodata) THEN
            h = dem % mat (ii+1,jj)
          ELSE
            h = dem % mat (ii,jj)
          END IF
        ELSE  
          h = dem % mat (ii,jj)  
        END IF
        
        !south-east cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj+1,dem) ) THEN
          IF ( dem % mat (ii+1,jj+1) /= dem % nodata) THEN
            i = dem % mat (ii+1,jj+1)
          ELSE
            i = dem % mat (ii,jj)
          END IF
        ELSE  
          i = dem % mat (ii,jj)  
        END IF
        
        !compute the rate of change in the x direction
        dZdX = ((c + 2.*f + i) - (a + 2.*d + g)) / 8.
        
        !compute the rate of change in the y direction
        
        dZdY = ((g + 2.*h + i) - (a + 2.*b + c))  / 8.
        
        !compute aspect
        IF (dZdX == 0. .AND. dZdY == 0.) THEN
            aspect % mat (ii,jj)  = -1.
            CYCLE
        ELSE
            aspectIJ  =  ATAN2 (dZdY, -dZdX)
        END IF
        
        
        !convert aspect value to compass direction values (0–2Pi radians)
        
        IF (aspectIJ < 0) THEN

            aspect % mat (ii,jj)  = pi/2. - aspectIJ

        ELSE IF (aspectIJ > pi /2.) THEN 

           aspect % mat (ii,jj) = 2.*pi - aspectIJ + pi/2.
           
        ELSE

            aspect % mat (ii,jj) = pi/2. - aspectIJ 
           
        END IF
       
    END IF
  END DO jdim
END DO idim 

RETURN

END SUBROUTINE DeriveAspect