0

我尝试根据以下代码进行 Cokrig，一切似乎都可以拟合 CoV，当我尝试生成插值时出现错误。这是我习惯于我的案例的代码：

library(terra)
library(gstat)
library(sp)
library(rgeos)
library(maptools)
library(rasterize)
library(raster)
library(lattice)
library(geoR)

setwd("C:/Users/c/Downloads/Documents")
Radata <- read.csv("raindata.csv", TRUE, ",")
ubngrid <- read.csv("UBNgridded.csv", TRUE, ",")
coordinates(Radata) = ~ X + Y
coordinates(ubngrid) = ~ X + Y

rast <- raster(ncol = 42, nrow = 49)
extent(rast) <- extent(ubngrid)
rasterize(ubngrid, rast, fun = mean)
class(ubngrid)
proj4string(ubngrid) <-CRS("+proj=utm +zone=37 +datum=WGS84")

gCoK <- gstat(NULL, 'Arsenic',Ar~1, Radata)
gCoK <- gstat(gCoK, 'Rainfall', obstrans~1, Radata)

coV <- variogram(gCoK)
plot(coV, type='b', main='Co-variogram')
coV.fit <- fit.lmc(coV, gCoK, vgm(model='Sph', range=1000))
coV.fit

plot(coV, coV.fit, main='Fitted Co-variogram')
coK <- interpolate(ubngrid,coV.fit,debug.level=0)
plot(coK)

错误是：

（函数（类，fdef，mtable）中的错误：无法找到签名“空间点”的函数“插值”的继承方法</p>

原始数据和网格的数据结构如下：

Raw data (Radata):
structure(list(X = c(292722.66, 292722.66, 347471.84, 331697.62, 
                     347308.09, 353000.68, 472622.24, 461595.36, 467048.61, 373573.47, 
                     329853.95, 221532.11, 330428.43, 319514.95, 314044.23, 324985.55, 
                     341345.11, 202653.33, 214451.07, 368598.76, 352241.31, 241704.19, 
                     258066.2, 236060.39, 357179.46, 5e+05, 170390.47, 111692.65, 
                     161189.76, 390750.98, 314386.44, 379824.88, 123283.36, 93047.35, 
                     280376.32, 291732.42, 445359.1, 280634.02, 192294.72, 236060.39
), Y = c(1266419.32, 1266419.32, 1304819.32, 1404460.8, 1271637.56, 
         1321383.37, 1133064.21, 1061207.17, 1022505.54, 989521.73, 1050511.46, 
         1239297.88, 1172173.82, 1177760.59, 1177790.16, 1177731.91, 1166591.08, 
         929542.94, 1040134.26, 1138827.91, 1155482.53, 1006746.05, 984515.35, 
         984652.09, 1017222.82, 1000380.05, 1029415.07, 1196121.75, 1217729.6, 
         1216161.36, 1238628.64, 1216199.57, 1251375.28, 996909.86, 1039701.79, 
         1106012.3, 1204967.51, 1083951.79, 1018166.55, 984652.09), obstrans = c(3.194, 
                                                                                 3.194, 3.164, 3.068, 3.102, 3.164, 2.948, 2.948, 3.076, 3.011, 
                                                                                 3.011, 2.987, 3.132, 3.132, 3.132, 3.132, 3.132, 3.316, 3.316, 
                                                                                 3.132, 3.132, 3.316, 3.316, 3.316, 3.143, 3.031, 3.052, 3.009, 
                                                                                 3.124, 3.124, 3.124, 3.124, 3.124, 3.157, 3.157, 3.142, 3.142, 
                                                                                 3.295, 3.316, 3.316), Ar = c(3.2, 5, 4.1, 0.5, 2.9, 0.2, 0.2, 
                                                                                                              1.2, 0.9, 2.5, 3.1, 0.7, 0.2, 0.3, 2.9, 2.1, 0.9, 2.6, 3, 1.7, 
                                                                                                              3.8, 7, 3.6, 2.9, 0.9, 2.5, 3.1, 0.7, 0.2, 0.3, 2.9, 2.1, 0.9, 
                                                                                                              2.6, 3, 1.7, 3.8, 7, 3.6, 2.9)), class = "data.frame", row.names = c(NA, 
                                                                                                                                                                                   -40L))

UBN grid data:
Grid profile summary: RasterLayer dimensions : 49, 42, 2058 (nrow, ncol, ncell) resolution : 9761.905, 9795.918 (x, y) extent : 93047, 503047, 929543, 1409543 (xmin, xmax, ymin, ymax) crs : +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs source : memory names : layer values : 1, 2058 (min, max)
4

1 回答 1

0

这是来自 的联合克里金示例?terra::interpolate。见?gstat和?raster::interpolate替代品

library(terra)
library(gstat)
meuse <- readRDS(system.file("ex/meuse.rds", package="terra"))
r <- rast(system.file("ex/meuse.tif", package="terra"))
gCoK <- gstat(NULL, 'log.zinc', log(zinc)~1, meuse, locations=~x+y)
gCoK <- gstat(gCoK, 'elev', elev~1, meuse, locations=~x+y)
gCoK <- gstat(gCoK, 'cadmium', cadmium~1, meuse, locations=~x+y)
gCoK <- gstat(gCoK, 'copper', copper~1, meuse, locations=~x+y)
coV <- variogram(gCoK)
plot(coV, type='b', main='Co-variogram')
coV.fit <- fit.lmc(coV, gCoK, vgm(model='Sph', range=1000))
coV.fit
plot(coV, coV.fit, main='Fitted Co-variogram')
coK <- interpolate(r, coV.fit, debug.level=0)
plot(coK)
于 2021-07-05T16:04:25.050 回答