DeriveSlope Subroutine

public subroutine DeriveSlope(dem, slope)

For each cell, the routine calculates the maximum rate of change in value from that cell to its neighbors. Basically, the maximum change in elevation over the distance between the cell and its eight neighbors identifies the steepest downhill descent from the cell. Conceptually, the routine fits a plane to the z-values of a 3 x 3 cell neighborhood around the processing or center cell. The slope value of this plane is calculated using the average maximum technique. The direction the plane faces is the aspect for the processing cell. If there is a cell location in the neighborhood with a NoData z-value, the z-value of the center cell will be assigned to the location.

References:

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

Arguments

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

Variables

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

length
real(kind=float), public :: a

elevation of N-W cell
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 DeriveSlope &
!
(dem, slope)


IMPLICIT NONE

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

!arguments with intent out
TYPE (grid_real), INTENT (out) :: slope ![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) :: 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
REAL (KIND = float) :: L !!length

!------------------------------end of declaration -----------------------------

!allocate slope grid
CALL NewGrid (slope,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
    
        !set cell size
        L = CellArea (dem,ii,jj) ** 0.5
        
        !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. * L)
        
        !compute the rate of change in the y direction
        
        dZdY = ((g + 2.*h + i) - (a + 2.*b + c))  / (8. * L)
        
        !compute slope
        slope % mat (ii,jj) = ATAN ( (dZdX ** 2. + dZdY ** 2.) ** 0.5 )
        
    
    END IF
  END DO jdim
END DO idim 


RETURN 
END SUBROUTINE DeriveSlope