继续之前的帖子，下述代码主要是针对地形分析。

1. 地形数据读取

esdf<-readxl::read_excel("data2.xlsx", sheet=1,skip=1)

2.地形因子生成im形式

##选出海拔的数据
 elev=esdf$meanelev 
 
##将海拔数据转换为30行50列的矩阵
 e=matrix(data=elev,nrow=30,ncol=50)
 
##x，y是每个样方中心的坐标
 x=seq(from=10,to=990,by=20)
 y=seq(from=10,to=590,by=20) 
 
##转换为im形式
 meanelev=im(e,x,y)         
 plot(meanelev)

3.bermantest函数用法--不太清楚是否是bermantest？

##函数1--as.list.vector()
## 在第二步S.cal函数中将list转化为向量，便于计算均值。
 
as.list.vector=function(x){ 
  ver=c()   
  for(i in 1:length(x)) ver[i]=x[[i]] 
  ver 
}

##函数2--cal.p.list()计算p的零模型p.list
##p代表物种实际分布的格局，为ppp格式
 cal.p.list=function(p){    
   ##最后p.list中包含1499个零模型和一个真实格局，均为ppp格式
   
   #这时p.list一个空列表，用来存放下一步生成的零模型
   p.list=list() 
   
   #行依次移动1到50个样方
 
 for(i in 1:50){
     #列依次移动1到30个样方，一共生成1500个零模型
     for(j in 1:30){          
       p.list[[(i-1)*30+j]]=ppp(tor.shift(as.points(p$x,p$y),xsh=20*i,
                                          ysh=20*j)[ ,1],
                                tor.shift(as.points(p$x,p$y),xsh=20*i,
                                          ysh=20*j)[ ,2],
                                c(0,1000),c(0,600))       
       
     }     
   } 
   return(p.list)
 }
 
##函数3--f()计算点（x,y）处的环境变量v值，
##其中xy是点的笛卡尔坐标；c.im是要测试的地形因子 
 f=function(xy,c.im){ 
   x=xy[1]#x坐标值
   y=xy[2]#y坐标轴
   rol.cal=function(y){
     if(floor(y/20)==30){floor(y/20)}else{floor(y/20)+1}}
   col.cal=function(x){
     if(floor(x/20)==50){floor(x/20)}else{floor(x/20)+1}}
   rol=rol.cal(y)
   col=col.cal(x)
   re=c.im[rol,col] ##用样方序号确定环境变量的值
   re
 }
##函数4--S.cal()计算点格局p的S值；
## cov是要测试的地形因子，对应于f函数中的c.im
 S.cal=function(p,cov){  
   xy.list=list()
   for (i in 1:npoints(p)){
     xy.list[[i]]=c(p[i]$x,p[i]$y)
   }

   v=lapply(xy.list,f,c.im=cov)
   S=mean(as.list.vector(v))
   S
 }
 
## 分析程序过程：
 sp.list.res<-cal.p.list(sp)

 ##Ssim.list是一个列表，包含1500个元素，其中有1499个Ssim和一个Sobs。
 Ssim.list=lapply(sp.list.res,S.cal,cov=meanelev) # very slow!  
 Ssim.ver=as.list.vector(Ssim.list)   
   
 ##计算真实格局的S值在这1500个中的排位
 n=rank(Ssim.ver)[1500] 
   
 if(n>1500*0.025&n<1500*(1-0.025)){"n"}else{ 
     if(n<=1500*0.025){"-"}else{"+"}
 }

运行结果：

## [1] "n"

最后的if结构表示：
   如果真实点格局的S值在这1500个值中排名在前0.025，返回值为“+”，在后0.025返回值为“-”，在0.025-0.975之间返回值为“n”。

4.bermantest函数用法--自带例子

# Berman's data

 X <- copper$SouthPoints
 L <- copper$SouthLines
 D <- distmap(L, eps=1)
# test of CSR
 berman.test(X, D)

运行结果：

 ## 
 ##  Berman Z1 test of CSR in two dimensions
 ## 
 ## data:  covariate 'D' evaluated at points of 'X'
 ## Z1 = 0.46962, p-value = 0.6386
 ## alternative hypothesis: two-sided