R语言：基于Spatstat包的单变量和双变量平面Ripley‘s K函数分析

Posted 2022-01-29 朱鷺子

因为网上查找得到的R语言K函数分析的资料很少，Spatstat包的介绍也很少，所以抛砖引玉。本文简单介绍了使用R语言spatstat包的Kest，Kcross(Lcross)和Kinhom方法进行了单变量和双变量平面Ripley's K函数分析的方法及代码，希望有所帮助。

首先需要将经纬度坐标转换为笛卡尔坐标系下的坐标，不然Ripley's K函数计算的空间尺度半径是不对的。

原始数据OriginalData的格式如下：

OriginalData.csv

OBJECTID	longitude	latitude
23	30.93427	121.5829
590	30.88591	121.5676
782	30.88255	121.5452
1085	30.93471	121.5365

 我们所需要的数据实际只有经纬度两列。读入OriginalData.csv后，可以通过SoDA包的geoXY方法进行转换，将结果保存在point.csv中。代码如下：

#将经纬度坐标转换为笛卡尔坐标
library("SoDA")
pos <- read.csv("OriginalData.csv", header = T)
geoxy <- geoXY(pos$longitude, pos$latitude)
write.csv(geoxy,file = "point.csv")

生成的结果文件point.csv如下：

point.csv

	X	Y
1	63623.28	23879.09
2	62155.98	18517.41
3	60007.67	18145
4	59173.4	23927.32

接下来先进行单变量平面Ripley's K函数分析，读入刚刚生成的笛卡尔坐标文件。可以看到Kest方法的关键是生成ppp对象，注意此处的范围window中使用了xrange和yrange，缺点是必须知道数据的上下限。我是在Excel中通过max()和min()取到数据的极值之后才来做的。也可以直接使用range函数，具体代码会补充在下面。

#进行单变量K函数分析：Spatstat包Kest方法
library(spatstat)
data1 <- read.csv("point.csv",header = T) #笛卡尔坐标系下的点坐标
df1 <- data1[,2:3]
sp <- ppp(df1$X, df1$Y, 
          window = owin(xrange=c(0,28770.79991),yrange=c(0,101406.9336)))
k <-Kest(sp)
plot(k,xlab="Scale/m",ylab="K(r)",main="单变量K函数分析")

直接使用range函数，就不需要事先获取数据的上下限，生成ppp对象的代码如下：

#使用range确定范围
point <- ppp(df1$X,df1$Y,
             range(df1$X),range(df1$Y))

最终可以得到plot绘制的K函数图像分析结果如下：

 摘录帮助文档中对Kest方法输出Value的描述：

r    the vector of values of the argument r at which the function K has been estimated

theo    the theoretical value K(r) = pi * r^2 for a stationary Poisson process

together with columns named "border", "bord.modif", "iso" and/or "trans", according to the selected edge corrections. These columns contain estimates of the function K(r) obtained by the edge corrections named.

接下来进行双变量的K函数计算。我手头的数据是两个点位数据，计算第一个点与第二个点的相关K函数。事先进行数据处理，将两类点位放在同一个文件position.csv中，用type进行标记。原始数据如下：

position.csv

ID	type	X	Y
1	station	11698.39	40730.17
2	station	12258.3	41581.32
…	…	…	…
479	point	63623.28	23879.09
480	point	62155.98	18517.41

 读入position.csv，此处的关键依然是生成ppp图像，注意需要用marks标记不同类的点数据。使用spatstat包的Kcross方法即可简单地完成双变量K函数计算，结果保存在Kab。

也可以使用Lcross进行L函数计算，结果保存在Lab。

#双变量K函数分析
data3 <- read.csv("position.csv")
df4 <- data3[,1:4]
df4$type <- as.factor(df4$type)
position <- ppp(df4$X,df4$Y,
              range(df4$X),range(df4$Y),marks = df4$type)
Kab <- Kcross(position, "station", "point",correction=c("Ripley","border"))
Lab <- Lcross(position, "station", "point",correction=c("Ripley","border"))
plot(Kab)

具体选择哪几类数据可以利用dplyr包的filter函数，比如选择type为point和station的数据时，代码如下：

positionsab <- dplyr::filter(df, type %in% c("point", "station"))

最终可以得到双变量平面K函数的图像：

 也可以利用Kinhom()进行双变量点的平面kinhom函数计算。

#双变量平面Kinhom函数分析
library("spatstat")
data <- read.csv("point.csv",header = T)
df <- data[,1:3]
point<-ppp(df$X,df$Y,
           range(df$X),range(df3$Y),marks = df$type)
#x和y坐标进行多项式拟合空间趋势函数
fit <- ppm(point, ~ polynom(x,y,2), Poisson())
lambda <- predict(fit, locations=point, type="trend")
Ki <- Kinhom(point, lambda)
plot(Ki,main="双变量平面Kinhom函数分析")

绘制图像结果：

以上是计算K函数的方法。如果要进行K函数的分析，需要通过envelope计算CSR随机分布模型，并比较空模型模拟的估计的上下限，如果在上限以上则表明是集聚分布，在下限以下则表明是均匀分布，在上下限之间则表明是随机分布。代码如下：

#读入数据
library(spatstat)
data <- read.csv("position.csv",header = T) #笛卡尔坐标系下的点坐标
df <- data1[,1:3]
point <-ppp(df$X,df$Y,
            range(df$X),range(df3$Y),marks = df$type)


#单变量K函数分析：Kest
pointsub <-dplyr::filter(df, type %in% c("point")) #筛选出类型为type的单变量点数据
pointsub$type <- as.factor(pointsub$type)
envelope1 <- envelope(pointsub,fun=Kest,nsim=100) #nsim表示迭代次数

plot(envelope1,xlab="Scale/m",ylab="K(r)",main="单变量点Ripley's K函数分析")


#双变量K函数分析：Kcross
envelope2 <- envelope(point,fun=Kcross,nsim=100)

plot(envelope2,xlab="Scale/m",ylab="K(r)",main="双变量点Ripley's K函数分析")


#双变量Kinhom函数分析
fit <- ppm(point, ~ polynom(x,y,2), Poisson()) #x和y坐标进行多项式拟合空间趋势函数
lambda <- predict(fit, locations=point, type="trend")
envelope3 <- envelope(fit,fun=Kinhom,nsim=100)

plot(envelope2,main="双变量Kinhom函数分析")

绘制的函数图像如下，hi和lo分别是CSR模型估计的上限和下限。

 

 最后，感谢林元震老师分享的生态R包spatstat的部分用法一文，获益匪浅！本文的代码也是基于林老师的工作基础。感谢大家阅读！