在 R 中进行克里金法之前,我应该绘制哪种变差函数?
What KInd of variogram should I plot before doing Kriging in R?
我有一个名为 seoul3112 的 csv 文件,包含 PM10 浓度。 please, download.。我尝试绘制样本变差函数并在其上拟合模型。
library(sp)
library(gstat)
library(rgdal)
seoul3112<-read.csv("seoul3112.csv",row.name=1)
seoul3112<-na.omit(seoul3112)
#assign a CRS and reproject
coordinates(seoul3112)=~LON+LAT
proj4string(seoul3112) = "+proj=longlat +datum=WGS84"
seoul3112<-spTransform(seoul3112, CRS("+proj=utm +north +zone=52 +datum=WGS84"))
#绘制半变异函数
g<-gstat(id="PM10",formula=PM10~LON+LAT, data=seoul3112)
seoul3112.var<-variogram(g, cutoff=70000, width=6000)
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",40000,400),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")
写完这些代码后,我得到了这样的半变异函数。
由于我在地质统计学方面还很陌生,所以我很困惑上面的变差函数是否适合我的数据集。因为,在典型的变异函数中,半方差值在基台处变为水平。但是这个变异函数向上!我应该在我的代码中做一些更正吗?
另一件事是,实际上我的最终目标是对我的数据集 (seoul3112) 进行克里金插值。我不明白,为了进行克里金法,这个样本变差函数就足够了,或者我应该绘制方向变差函数或任何其他东西?
谁能详细解释一下吗?
如果你查看你的拟合模型,它会有一个 sill
参数(即 nugget
+ psill
),但它超出了你的样本范围。
sill = sum(model.3112$psill)
您可能没有足够远的点对距离到达 range
到 sill
。我认为这不是问题,只要您使用此半变异函数在数据的空间域 + 60000 米(您有数据支持的距离)内进行预测。我采取这个立场是因为最重要的半变异函数拟合部分是开始,并且在数据范围内拟合很好。
要检查的是有多少点对 (np
) 支持您绘制的 bins
(点对越多越好)。
另一个建议是使用水平图寻找各向异性
seoul3112.var_map<-variogram(g, cutoff=70000, width=6000, map=TRUE)
plot(seoul3112.var_map, col.regions=terrain.colors(20))
或者通过设置 alpha
并绘制模型拟合来查看多个方向的拟合。您只需检查 180 度,因为它是对称的(您可以在下图中看到这一点,其中 0 和 180 是相同的)。
seoul3112.var_angles<-variogram(g, cutoff=70000, width=6000, alpha=seq(0,180,30))
plot(seoul3112.var_angles, model=model.3112, pch=16, ylim=c(0,3000))
如果在方向上存在差异,您可以使用 fit.variogram
为各向异性的长轴和短轴建模,并为 vgm()
模型定义了 anis
。
示例:
vgm(700,"Gau",40000,400, anis=c(someangle, someratio))
我有一个名为 seoul3112 的 csv 文件,包含 PM10 浓度。 please, download.。我尝试绘制样本变差函数并在其上拟合模型。
library(sp)
library(gstat)
library(rgdal)
seoul3112<-read.csv("seoul3112.csv",row.name=1)
seoul3112<-na.omit(seoul3112)
#assign a CRS and reproject
coordinates(seoul3112)=~LON+LAT
proj4string(seoul3112) = "+proj=longlat +datum=WGS84"
seoul3112<-spTransform(seoul3112, CRS("+proj=utm +north +zone=52 +datum=WGS84"))
#绘制半变异函数
g<-gstat(id="PM10",formula=PM10~LON+LAT, data=seoul3112)
seoul3112.var<-variogram(g, cutoff=70000, width=6000)
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",40000,400),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")
写完这些代码后,我得到了这样的半变异函数。
由于我在地质统计学方面还很陌生,所以我很困惑上面的变差函数是否适合我的数据集。因为,在典型的变异函数中,半方差值在基台处变为水平。但是这个变异函数向上!我应该在我的代码中做一些更正吗?
另一件事是,实际上我的最终目标是对我的数据集 (seoul3112) 进行克里金插值。我不明白,为了进行克里金法,这个样本变差函数就足够了,或者我应该绘制方向变差函数或任何其他东西? 谁能详细解释一下吗?
如果你查看你的拟合模型,它会有一个 sill
参数(即 nugget
+ psill
),但它超出了你的样本范围。
sill = sum(model.3112$psill)
您可能没有足够远的点对距离到达 range
到 sill
。我认为这不是问题,只要您使用此半变异函数在数据的空间域 + 60000 米(您有数据支持的距离)内进行预测。我采取这个立场是因为最重要的半变异函数拟合部分是开始,并且在数据范围内拟合很好。
要检查的是有多少点对 (np
) 支持您绘制的 bins
(点对越多越好)。
另一个建议是使用水平图寻找各向异性
seoul3112.var_map<-variogram(g, cutoff=70000, width=6000, map=TRUE)
plot(seoul3112.var_map, col.regions=terrain.colors(20))
或者通过设置 alpha
并绘制模型拟合来查看多个方向的拟合。您只需检查 180 度,因为它是对称的(您可以在下图中看到这一点,其中 0 和 180 是相同的)。
seoul3112.var_angles<-variogram(g, cutoff=70000, width=6000, alpha=seq(0,180,30))
plot(seoul3112.var_angles, model=model.3112, pch=16, ylim=c(0,3000))
如果在方向上存在差异,您可以使用 fit.variogram
为各向异性的长轴和短轴建模,并为 vgm()
模型定义了 anis
。
示例:
vgm(700,"Gau",40000,400, anis=c(someangle, someratio))