R中如何对spatialEco包的O-ring统计量曲线积分求曲线下面积

问题背景 #

在RStudio中计算空间点格局的O-ring统计量时，spatialEco包的返回结果不像spatstat包中Ripley's K函数那样是标准化的fv（函数值）对象，因此无法直接通过通用方法调取原始计算值，需要基于输出结构实现曲线下方面积估算。

当前实现代码 #

library('spatstat')
library('ggplot2')
library('spatialEco')
library(dplyr) # 原代码使用slice函数做子集筛选，需加载该包才能正常运行

set.seed(seed=24)

radiusCluster<-100
lambdaParent<-.02
lambdaDaughter<-30
hosts<-1000
randmod<-0
dim<-2000


numbparents<-rpois(1,lambdaParent*dim)

xxParent<-runif(numbparents,0+radiusCluster,dim-radiusCluster)
yyParent<-runif(numbparents,0+radiusCluster,dim-radiusCluster)

numbdaughter<-rpois(numbparents,(lambdaDaughter))
sumdaughter<-sum(numbdaughter)


thetaLandscape<-2*pi*runif(sumdaughter)

rho<-radiusCluster*sqrt(runif(sumdaughter))


xx0=rho*cos(thetaLandscape)
yy0=rho*sin(thetaLandscape)


xx<-rep(xxParent,numbdaughter)
yy<-rep(yyParent,numbdaughter)

xx<-xx+xx0

yy<-yy+yy0
cds<-data.frame(xx,yy)
is_outlier<-function(x){
  x > dim| x < 0
}
cds<-cds[!(is_outlier(cds$xx)|is_outlier(cds$yy)),]
while (nrow(cds)<hosts){
 dif<-hosts-nrow(cds)
  extraparentxx<-sample(xxParent,dif,replace = TRUE)
  extraparentyy<-sample(yyParent,dif,replace = TRUE)
  extrathetaLandscape<-2*pi*runif(dif)
  extrarho<-radiusCluster*sqrt(runif(dif))
  newextracoodsxx<-extrarho*cos(extrathetaLandscape)
  newextracoodsyy<-extrarho*sin(extrathetaLandscape)
  extraxx<-extraparentxx+newextracoodsxx
  extrayy<-extraparentyy+newextracoodsyy
  cdsextra<-data.frame(xx=extraxx,yy=extrayy)
  cds<-rbind(cds,cdsextra)
}


sampleselect<-sample(1:nrow(cds),hosts,replace=F)
cds<-cds%>%slice(sampleselect)

randfunction<-function(x){
  x<-runif(length(x),0,dim)
}
randselect<-sample(1:nrow(cds),floor(hosts*randmod),replace=F)
cds[randselect,]<-apply(cds[randselect,],1,randfunction)

landscape<-ppp(x=cds$xx,y=cds$yy,window=owin(xrange=c(0,dim),yrange=c(0,dim)))
ggplot(data.frame(landscape))+geom_point(aes(x=x,y=y))+coord_equal()+theme_minimal()

ostat<-o.ring(landscape,inhomogenous=FALSE)

运行上述代码后生成的O-ring统计图如下：

曲线下面积计算实现 #

1. 提取O-ring统计量原始值 #

spatialEco包的o.ring函数返回的结果为结构化列表，不需要从渲染的图中逆向提取数据，直接通过列表索引即可获取全部分离的计算值：

• ostat$r：存储所有参与计算的半径距离序列（对应曲线的横轴值）
• ostat$v：存储每个半径对应的O-ring统计量值（对应曲线的纵轴值）
  若索引报错，可执行str(ostat)查看返回对象的实际字段名，根据输出结果调整索引即可

2. 数值积分计算面积 #

提取到横轴、纵轴序列后，采用梯形数值积分法计算曲线下面积即可：

• 可加载caTools包，使用内置的trapz函数直接计算：

# 首次运行需先安装包
install.packages("caTools")
library(caTools)

# 计算全段曲线下面积
oring_auc <- trapz(x = ostat$r, y = ostat$v)

• 若不想额外安装包，可手动实现梯形积分逻辑：遍历所有相邻的点对，累加每个相邻点构成的梯形面积，即可得到总积分值。
• 若仅需要计算某一特定半径区间内的面积，可先筛选出落在目标半径区间内的r和对应O-ring值，再执行积分操作即可。

内容的提问来源于stack exchange，提问作者OpenSauce