|||
以下函数, 用来计算森林监测样地每个小样方的坡度和坡向以及凹凸度。 这里用到的是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)
}
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 11:23
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社