yzhlinscau的个人博客分享 http://blog.sciencenet.cn/u/yzhlinscau

博文

生态R包spatstat的部分用法(三)--地形分析

已有 4564 次阅读 2019-1-7 11:20 |个人分类:R|系统分类:科研笔记| 空间生态学, 地形分析

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


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)

blob.png


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


https://blog.sciencenet.cn/blog-1114360-1155752.html

上一篇:生态R包spatstat的部分用法(二)
下一篇:ASReml-R与动植物遗传评估实训会正式通知(2019年4月7日-10日,广州)
收藏 IP: 125.88.24.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 19:46

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部