样地数据从海拔计算坡度、坡向和凹凸度

以下函数， 用来计算森林监测样地每个小样方的坡度和坡向以及凹凸度。 这里用到的是Arcgis中的计算方法。
各函数都假设，原点位于西南方向， x向东增加， y向北增加。 函数都已经经过演算， 结果与示例结果相同。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
## 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)
}