R语言：克里金插值

 

基于空间自相关，R语言克里金插值

library(gstat)

Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string

library(sp)
data(meuse)
head(meuse)

       x      y cadmium copper lead zinc  elev       dist   om ffreq soil lime landuse
1 181072 333611    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1      Ah
2 181025 333558     8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1      Ah
3 181165 333537     6.5     68  199  640 7.800 0.10302900 13.0     1    1    1      Ah
4 181298 333484     2.6     81  116  257 7.655 0.19009400  8.0     1    2    0      Ga
5 181307 333330     2.8     48  117  269 7.480 0.27709000  8.7     1    2    0      Ah
6 181390 333260     3.0     61  137  281 7.791 0.36406700  7.8     1    2    0      Ga
  dist.m
1     50
2     30
3    150
4    270
5    380
6    470

数据框转空间数据

coordinates(meuse) <- c("x","y")#定义坐标 
class(meuse)

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

spplot(meuse,"cadmium")

In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

v

    np       dist     gamma dir.hor dir.ver   id
1   57   79.29244 0.6650872       0       0 var1
2  299  163.97367 0.8584648       0       0 var1
3  419  267.36483 1.0064382       0       0 var1
4  457  372.73542 1.1567136       0       0 var1
5  547  478.47670 1.3064732       0       0 var1
6  533  585.34058 1.5135658       0       0 var1
7  574  693.14526 1.6040086       0       0 var1
8  564  796.18365 1.7096998       0       0 var1
9  589  903.14650 1.7706890       0       0 var1
10 543 1011.29177 1.9875659       0       0 var1
11 500 1117.86235 1.8259154       0       0 var1
12 477 1221.32810 1.8852099       0       0 var1
13 452 1329.16407 1.9145967       0       0 var1
14 457 1437.25620 1.8505336       0       0 var1
15 415 1543.20248 1.8523791       0       0 var1

np点对数量，dist距离，gamma两点协方差

半变异函数画图

plot(v,plot.number = T)

设置半变异计算时点对距离

v <- variogram(log(meuse$cadmium) ~ 1,meuse, width = 100)
plot(v, plot.number = T)

半变异函数值越大协方差越小。Nugget表示不确定性，Sill最大gamma值，Range半变异gamma不再增加的点。 PSill= sill - nugget. 球形曲线、高斯曲线

经验半变异函数对象

m1 <- vgm(psill = 1.5,
          model = "Sph",#Sph球形曲线
          range = 1400,
          nugget = 0.5)

Warning messages:
1: In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string
2: In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

m2 <- fit.variogram(v,m1)#拟合半变异函数"

plot(v,model = m2)

拟合计算出的模型参数

m2

  model     psill    range
1   Nug 0.5855785    0.000
2   Sph 1.3253037 1224.936

开始空间插值

数据准备

data(meuse.grid)
data(meuse)
head(meuse.grid)

       x      y part.a part.b      dist soil ffreq
1 181180 333740      1      0 0.0000000    1     1
2 181140 333700      1      0 0.0000000    1     1
3 181180 333700      1      0 0.0122243    1     1
4 181220 333700      1      0 0.0434678    1     1
5 181100 333660      1      0 0.0000000    1     1
6 181140 333660      1      0 0.0122243    1     1

普通克里金插值

kr <- krige(log(meuse$cadmium) ~ 1,
            loc = ~ x+y,            #位置信息
            data = meuse,           
            newdata = meuse.grid,   
            model = m2)

[using ordinary kriging]

head(kr)

       x      y var1.pred  var1.var
1 181180 333740  1.831241 1.1606001
2 181140 333700  1.905750 1.0508440
3 181180 333700  1.834357 1.0788180
4 181220 333700  1.760312 1.1129144
5 181100 333660  1.985888 0.9396745
6 181140 333660  1.909244 0.9662369

var1.var 估计值的方差，越大越不可靠。

显示插值结果图

In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

levelplot(var1.pred ~ x + y,kr )

地理探测器，基于异质性，R语言插值

http://www.geodetector.org/

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQojIOWfuuS6juepuumXtOiHquebuOWFs++8jFLor63oqIDlhYvph4zph5Hmj5LlgLwNCg0KYGBge3J9DQpsaWJyYXJ5KGdzdGF0KQ0KbGlicmFyeShzcCkNCmRhdGEobWV1c2UpDQpoZWFkKG1ldXNlKQ0KDQpgYGANCiMjIOaVsOaNruahhui9rOepuumXtOaVsOaNrg0KYGBge3J9DQpjb29yZGluYXRlcyhtZXVzZSkgPC0gYygieCIsInkiKSPlrprkuYnlnZDmoIcgDQpjbGFzcyhtZXVzZSkNCmBgYA0KDQpgYGB7cn0NCnNwcGxvdChtZXVzZSwiY2FkbWl1bSIpDQpgYGANCmBgYHtyfQ0KdiA8LSB2YXJpb2dyYW0obG9nKG1ldXNlJGNhZG1pdW0pIH4gMSxtZXVzZSkj5Y2K5Y+Y5byC5Ye95pWwDQp2DQpgYGANCm5w54K55a+55pWw6YeP77yMZGlzdOi3neemu++8jGdhbW1h5Lik54K55Y2P5pa55beuDQoNCiMjIOWNiuWPmOW8guWHveaVsOeUu+Wbvg0KYGBge3J9DQpwbG90KHYscGxvdC5udW1iZXIgPSBUKQ0KYGBgDQojIyMg6K6+572u5Y2K5Y+Y5byC6K6h566X5pe254K55a+56Led56a7DQpgYGB7cn0NCnYgPC0gdmFyaW9ncmFtKGxvZyhtZXVzZSRjYWRtaXVtKSB+IDEsbWV1c2UsIHdpZHRoID0gMTAwKQ0KcGxvdCh2LCBwbG90Lm51bWJlciA9IFQpDQpgYGANCuWNiuWPmOW8guWHveaVsOWAvOi2iuWkp+WNj+aWueW3rui2iuWwj+OAgk51Z2dldOihqOekuuS4jeehruWumuaAp++8jFNpbGzmnIDlpKdnYW1tYeWAvO+8jFJhbmdl5Y2K5Y+Y5byCZ2FtbWHkuI3lho3lop7liqDnmoTngrnjgIINClBTaWxsPSBzaWxsIC0gbnVnZ2V0LiDnkIPlvaLmm7Lnur/jgIHpq5jmlq/mm7Lnur8NCg0KIyMg57uP6aqM5Y2K5Y+Y5byC5Ye95pWw5a+56LGhDQpgYGB7cn0NCm0xIDwtIHZnbShwc2lsbCA9IDEuNSwNCiAgICAgICAgICBtb2RlbCA9ICJTcGgiLCNTcGjnkIPlvaLmm7Lnur8NCiAgICAgICAgICByYW5nZSA9IDE0MDAsDQogICAgICAgICAgbnVnZ2V0ID0gMC41KQ0KbTIgPC0gZml0LnZhcmlvZ3JhbSh2LG0xKSPmi5/lkIjljYrlj5jlvILlh73mlbANCg0KYGBgDQpgYGB7cn0NCnBsb3Qodixtb2RlbCA9IG0yKQ0KYGBgDQoNCiMjIyDmi5/lkIjorqHnrpflh7rnmoTmqKHlnovlj4LmlbANCmBgYHtyfQ0KbTINCmBgYA0KIyDlvIDlp4vnqbrpl7Tmj5LlgLwNCg0KIyMg5pWw5o2u5YeG5aSHDQpgYGB7cn0NCmRhdGEobWV1c2UuZ3JpZCkNCmRhdGEobWV1c2UpDQpoZWFkKG1ldXNlLmdyaWQpDQpgYGANCiMjIOaZrumAmuWFi+mHjOmHkeaPkuWAvA0KYGBge3J9DQprciA8LSBrcmlnZShsb2cobWV1c2UkY2FkbWl1bSkgfiAxLA0KICAgICAgICAgICAgbG9jID0gfiB4K3ksICAgICAgICAgICAgI+S9jee9ruS/oeaBrw0KICAgICAgICAgICAgZGF0YSA9IG1ldXNlLCAgICAgICAgICAgDQogICAgICAgICAgICBuZXdkYXRhID0gbWV1c2UuZ3JpZCwgICANCiAgICAgICAgICAgIG1vZGVsID0gbTIpDQpoZWFkKGtyKQ0KYGBgDQp2YXIxLnZhciDkvLDorqHlgLznmoTmlrnlt67vvIzotorlpKfotorkuI3lj6/pnaDjgIINCg0KIyMg5pi+56S65o+S5YC857uT5p6c5Zu+DQpgYGB7cn0NCmxpYnJhcnkobGF0dGljZSkj55S75Zu+5YyFDQpsZXZlbHBsb3QodmFyMS5wcmVkIH4geCArIHksa3IgKQ0KYGBgDQojIOWcsOeQhuaOoua1i+WZqO+8jOWfuuS6juW8gui0qOaAp++8jFLor63oqIDmj5LlgLwNCg0KaHR0cDovL3d3dy5nZW9kZXRlY3Rvci5vcmcvDQoNCg0KDQo=