数据预处理

这次作业我用Edzer Pebesma,2003年开发的gstat包，去做变异函数，把老师给的图片转换为了数据,例如(x,y)=(3，4)表示从左下角开始，向右向上的第3列，第4行；删除了NA值，并运用coordinates函数建立了坐标系，将数据转换为空间数据.

Gräler B, Pebesma E, Heuvelink G (2016). “Spatio-Temporal Interpolation using gstat.” The R Journal, 8, 204-218. https://journal.r-project.org/archive/2016/RJ-2016-014/index.html.

library(gstat)
library(sp)
data<-read.csv("tata2.csv")
coordinates(data)=~x+y
head(data)
##   coordinates  z
## 1      (1, 1) 38
## 2      (1, 2) 36
## 3      (1, 3) 35
## 4      (1, 4) 37
## 5      (1, 5) 42
## 6      (1, 6) 40

空间数据作图

原来的数据

原来的数据

可视化

spplot(data,"z")

求变异函数值

np：the number of point pairs for this estimate

dist:the average distance of all point pairs considered for this estimate

gamma:the actual sample variogram estimate

(a<-variogram(z~1, data =data))
##    np     dist     gamma dir.hor dir.ver   id
## 1  72 1.000000  3.402778       0       0 var1
## 2  63 1.414214  5.500000       0       0 var1
## 3  60 2.000000  5.791667       0       0 var1
## 4 107 2.236068  6.672897       0       0 var1
## 5  43 2.828427  8.872093       0       0 var1
## 6  48 3.000000 10.104167       0       0 var1

半变异函数作图

拟合半变异函数

m1 <- vgm(psill = 0.5,
          model = "Sph",#Sph球形曲线
          range = 3,
          nugget = 0.5)
v <- variogram(z~ 1,coordinates(data),cutoff=10,width=1,data)
m2 <- fit.variogram(v,m1)
plot(v,model=m2,plot.number = T)

四个方向的变异函数

alpha：

direction in plane (x,y), in positive degrees clockwise from positive y (North): alpha=0 for direction North (increasing y), alpha=90 for direction East (increasing x); optional a vector of directions in (x,y)

  • alpha=0:南北
  • alpha=45:东北-西南
  • alpha=90:东西
  • alpah=135:西北-东南
x<-variogram(z~1,data =data,alpha=c(0,45,90,135))
x
##    np     dist     gamma dir.hor dir.ver   id
## 1  36 1.000000  5.347222       0       0 var1
## 2  27 2.000000  9.129630       0       0 var1
## 3  21 3.000000 17.547619       0       0 var1
## 4  31 1.414214  3.887097      45       0 var1
## 5  52 2.236068  4.644231      45       0 var1
## 6  22 2.828427  5.886364      45       0 var1
## 7  36 1.000000  1.458333      90       0 var1
## 8  33 2.000000  3.060606      90       0 var1
## 9  27 3.000000  4.314815      90       0 var1
## 10 32 1.414214  7.062500     135       0 var1
## 11 55 2.236068  8.590909     135       0 var1
## 12 21 2.828427 12.000000     135       0 var1
plot(x)