以下函数， 用来计算森林监测样地每个小样方的坡度和坡向以及凹凸度。 这里用到的是Arcgis中的计算方法。

各函数都假设，原点位于西南方向， x向东增加， y向北增加。 函数都已经经过演算， 结果与示例结果相同。

## Calculating theslopes for each grid in a forest dynamic plot

## Provide thematrix of altitudes

## Provide thecell_size, the altitudes and cell_size must be in the same unit.

## Algorithms couldbe found at ##

##http://webhelp.esri.com/arcgiSDEsktop/9.3/body.cfm?tocVisable=1&ID=-1&TopicName=How%20Slope%20works

## Value: a matrixrepresenting the degrees of each grid.

get_slope <-function(z.mat,cell_size =5){

    ###Extent the z.mat matrix using the value from the border, so that the slope forthe out most cells could be calculated.

    z.mat.colbind <-cbind(z.mat[,1], z.mat,z.mat[,ncol(z.mat)])

    z.mat.rbind <-rbind(c(z.mat[1,1], z.mat[1,], z.mat[1,ncol(z.mat)]),z.mat.colbind,

                         c(z.mat[nrow(z.mat),1], z.mat[nrow(z.mat),], z.mat[nrow(z.mat), ncol(z.mat)]))

    z.mat.ext <- z.mat.rbind

    nrz <-nrow(z.mat.ext)

    ncz <-ncol(z.mat.ext)

    slope.mat <-rep(NA, nrz*ncz)

    dim(slope.mat)<-c(nrz, ncz)  ##Create a Null matrix used in loop

    for(m in2:(nrz-1)){

       for(n in2:(ncz-1)){

           a = z.mat.ext[m-1,n-1]

           b = z.mat.ext[m-1,n]

           c=z.mat.ext[m-1,n+1]

           d = z.mat.ext[m,n-1]

           e = z.mat.ext[m,n]

           f = z.mat.ext[m,n+1]

           g = z.mat.ext[m+1,n-1]

           h = z.mat.ext[m+1,n]

           i = z.mat.ext[m+1,n+1]

           dz_dx =((c+2*f + i)-(a +2*d + g))/(8*cell_size)

           dz_dy =((g +2*h + i)-(a +2*b +c))/(8*cell_size)

           slope.mat[m,n]<-atan(sqrt( dz_dx^2+ dz_dy^2))*(180/pi)  

       }

   }

    return(slope.mat[c(-1, -nrow(slope.mat)),c(-1, -ncol(slope.mat))])

}

####

#### test.dat2 <-c(50, 30, 8, 45, 30, 10, 50, 30, 10)

#### test.dat.slope2<- matrix(test.dat2, nrow = 3, byrow = FALSE)

#### test.dat.slope2

####get_slope(test.dat.slope2, cell_size = 5)

###################################################################

## Calculating theaspect of each cell within a forest dynamic plot

## Adopted from ESRIhttp://webhelp.esri.com/arcgisdesktop/9.3/body.cfm?tocVisable=1&ID=-1&TopicName=how%20aspect%20works

## the z.mat is anumerical matrix representing the altitude for each cell

## the output is anumerical matrix representing the degree of aspect. Starting from the north,

## This functionassums that The direction of the column, from greater to smaller columns pointsto the north.

get_aspect <-function(z.mat){

    z.mat.colbind <-cbind(z.mat[,1], z.mat,z.mat[,ncol(z.mat)])

    z.mat.rbind <-rbind(c(z.mat[1,1], z.mat[1,], z.mat[1,ncol(z.mat)]),z.mat.colbind,

                         c(z.mat[nrow(z.mat),1], z.mat[nrow(z.mat),], z.mat[nrow(z.mat), ncol(z.mat)]))

    z.mat.ext <- z.mat.rbind

    nrz <-nrow(z.mat.ext)

    ncz <-ncol(z.mat.ext)

    aspect.mat <-rep(NA, nrz*ncz)

    dim(aspect.mat)<-c(nrz, ncz)

    for(m in2:(nrz-1)){

       for(n in2:(ncz-1)){

           a = z.mat.ext[m-1,n-1]

           b = z.mat.ext[m-1,n]

           c=z.mat.ext[m-1,n+1]

           d = z.mat.ext[m,n-1]

           e = z.mat.ext[m,n]

           f = z.mat.ext[m,n+1]

           g = z.mat.ext[m+1,n-1]

           h = z.mat.ext[m+1,n]

           i = z.mat.ext[m+1,n+1]

           dz_dx =((c+2*f + i)-(a +2*d + g))/8

           dz_dy =((g +2*h + i)-(a +2*b +c))/8

           aspect.mat[m, n]=(180/pi)* atan2(dz_dy, -dz_dx)

       }

   }

       format.aspect <-function(aspect){

          cell = aspect

          for(i in1:length(aspect)){

          if(aspect[i]<0){        

             cell[i]=90.0- aspect[i]        

             }        

          elseif(aspect[i]>90.0){        

             cell[i]=360.0- aspect[i]+90.0        

             }        

          else{        

             cell[i]=90.0- aspect[i]        

             }        

          }

          return(cell)        

       }

   aspect.mat <- aspect.mat[c(-1, -nrow(aspect.mat)),c(-1, -ncol(aspect.mat))]

   return(format.aspect(aspect.mat))

}

############################################################################

#### test.dat.aspect<- c(101, 101, 101, 92, 92, 91, 85, 85, 84)

####dim(test.dat.aspect) <- c(3, 3)

####image(test.dat.aspect)

#### test.dat.aspect

####get_aspect(test.dat.aspect)

########################Calculating  Convexity###################

get_convexity <-function(z.mat){

    z.mat.colbind <-cbind(rep(NA, nrow(z.mat)), z.mat, rep(NA, nrow(z.mat)))

    z.mat.rbind <-rbind(rep(NA, ncol(z.mat)+2), z.mat.colbind, rep(NA, ncol(z.mat)+2))

    z.mat.ext <- z.mat.rbind

    nrz <-nrow(z.mat.ext)

    ncz <-ncol(z.mat.ext)

    convexity.mat <-rep(NA, nrz*ncz)

    dim(convexity.mat)<-c(nrz, ncz)

    for(m in2:(nrz-1)){

       for(n in2:(ncz-1)){

           a = z.mat.ext[m-1,n-1]

           b = z.mat.ext[m-1,n]

           c=z.mat.ext[m-1,n+1]

           d = z.mat.ext[m,n-1]

           e = z.mat.ext[m,n]

           f = z.mat.ext[m,n+1]

           g = z.mat.ext[m+1,n-1]

           h = z.mat.ext[m+1,n]

           i = z.mat.ext[m+1,n+1]

           #### Ifna.....

           all.neighbour <-c(a, b, c, d, f, g, h, i)

           mean.neighbour <-mean(na.omit(all.neighbour))

           convexity.mat[m,n]<-(e - mean.neighbour)

       }

   }

   convexity.mat <- convexity.mat[c(-1, -nrow(convexity.mat)),c(-1, -ncol(convexity.mat))]

   return(convexity.mat)

}