基于空间自相关，R语言克里金插值
library(gstat)Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
EOF within quoted stringlibrary(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 localev 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 var1np点对数量，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 localem2 <- 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.9662369var1.var 估计值的方差，越大越不可靠。
显示插值结果图
In strsplit(code, "\n", fixed = TRUE) :
input string 1 is invalid in this localelevelplot(var1.pred ~ x + y,kr )地理探测器，基于异质性，R语言插值
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=
请发表评论