问题描述
我已尝试基于以下代码进行 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",")
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,debug.level=0)
plot(coK)
错误是：
（函数（类，fdef，mtable）中的错误： 无法为签名“SpatialPoints”找到函数“interpolate”的继承方法
原始数据和网格的数据结构如下：
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,2.948,3.076,3.011,2.987,3.132,3.316,3.143,3.031,3.052,3.009,3.124,3.157,3.142,3.295,3.316),Ar = c(3.2,5,4.1,0.5,2.9,0.2,1.2,0.9,2.5,3.1,0.7,0.3,2.1,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)
解决方法
这是一个来自 ?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,'cadmium',cadmium~1,'copper',copper~1,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,debug.level=0)
plot(coK)
如果您无法针对您的数据进行调整，请提供一个最小的、可重现的、独立的示例。我看到您在评论中包含了一些数据，但您应该编辑您的问题，以便我们可以复制您的代码并运行它（如上例所示）。另外，请删除识别问题和回答您的问题不需要的任何内容。