首页 文章

如何在R中使用gstat拟合半变异函数模型？

提问于
浏览 次
1

我有一个csv文件文件包含3月1日下午12点的大气PM10浓度数据 . please, download.我想在R中使用gstat包绘制一个半变异函数 . 我试着在R中编写这些代码 . 但是对于这些数据，我不能适应模型 .

library(sp)
    library(gstat)

    seoul3112<-read.csv("seoul3112.csv")
    seoul3112<-na.omit(seoul3112)

    g<-gstat(id="PM10",formula=PM10~LON+LAT,location=~LON+LAT,
             data=seoul3112)
    seoul3112.var<-variogram(g,width=0.04,cutoff=0.6)
    seoul3112.var
    plot(seoul3112.var, col="black", pch=16,cex=1.3,
         xlab="Distance",ylab="Semivariance",
         main="Omnidirectional Variogram for seoul 3112")

model.3112<- fit.variogram(seoul3112.var,vgm(700,"Gau",0.5,200), fit.method = 2)
    plot(seoul3112.var,model=model.3112, col="black", pch=16,cex=1.3,
         xlab="Distance",ylab="Semivariance",
         main="Omnidirectional Variogram for seoul 3112")

其实我是R和统计学的初学者 . 所以，即使是变差函，我也很无知 . 我有一些疑问：

a）当我将数据绘制为半变异函数时，它看起来不同于典型的半变异函数！为什么会这样？我应该对我的数据做任何其他改变吗？

b）如何使用这些数据拟合模型？我尝试了不同的模型，如“Sph”，“Exp”，但它们看起来像线性！为什么？

c）我怎样才能理解我应该在 vgm() 函数中使用的窗台，范围，金块的初始值？

d）我如何理解该模型是否适合数据？

e）对于使用克里金法，我应该绘制什么样的半变异函数？只有全方位半变异函数？或者我应该绘制方向半变异函数？

f）我怎样才能解释半变异函数？我的意思是从半变异函数的数据中我能理解的是什么？

提前致谢 .

r spatial kriging gstat covariogram

2 回答

  • 2

    我将为您的代码相关问题提供答案 . 其余的问题（d，e和f）更多地与理论相关 .

    首先，在您的评论中，当您更改 proj4string 时，距离单位应该在图上更改 . 他们做了吗？根据你的评论，听起来似乎没有发生 .

    a）除了玩 cutoff 距离之外，还要注意在半变异函数上支持每个 bin 的 np （点对） . 例如，使用您更新的 proj4string 信息，我尝试了 cutoff=80 和 width=80/10 （10个箱而不是15个）来查看半变异函数形状如何变化 . 在点对存在的情况下，从15到10个区域的减少不会改变，只会增加每个区域的距离 . 此外，这种方法不一定是您应该使用的方法，但它是如何更改二进制文件以获得更平滑的样本半变异函数的一个示例（但更平滑并不意味着更好） .

    Comparison of 10 Bins to 15 Bins

    b）使用您的代码， "Sph" 和 "Exp" 模型返回 Warning: singular model in variogram fit . 该警告表明没有足够的数据来拟合球形和指数经验模型的某些参数 . 有关每个经验公式及其参数的指导，请参见gstat user manual .

    c）例如， vgm() 函数可用于眼睛拟合样本半变异函数 . 如果您对如何使用样本数据绘制 vgm() 模型感到困惑，请尝试类似的方法

    eye_vgm = vgm(psill=1200,model="Gau",range=60,nugget=350)
plot(seoul3112.var,model=eye_vgm, col="black", pch=16,cex=1.3)

    您在 fit.variogram() 的调用中使用 vgm() ，因此只要您给 vgm() 的参数合理（例如基于样本数据）并且经验模型可以拟合参数， fit.variogram() 将根据 fit.method 找到拟合 .

    回复于
  • 0

    你的坐标是纬度和经度，但你不告知gstat它们是 . 因此，gstat将假设它可以计算这些数字的欧几里德距离，这是没有意义的 .

    建议是在使用包 sp 将您的点转换为 SpatialPointsDataFrame 之后学习如何使用gstat，然后学习如何投影数据以使欧几里德距离有意义 .

    回复于

相关问题