步骤：
1）载入原始的obs,discrete file # load the original obs.discrete filelibrary(atakrig)library(rgdal)## load demo data from rtop package ----if (!require("rtop", quietly = TRUE)) message("rtop library is required for demo data.")rpath <- system.file("extdata", package="rtop")observations <- readOGR(rpath, "observations")observations$obs <- observations$QSUMMER_OB/observations$AREASQKM## point-scale variogram ----obs.discrete <- discretizePolygon(observations, cellsize=1500, id="ID", value="obs")pointsv <- deconvPointVgm(obs.discrete, model="Exp", ngroup=12, rd=0.75, fig=TRUE)## cross validation ----pred.cv <- ataKriging.cv(obs.discrete, nfold=length(observations), pointsv)names(pred.cv)[6] <- "obs"summary(pred.cv[,c("obs","pred","var")])cor(pred.cv$obs, pred.cv$pred)# Pearson correlationmean(abs(pred.cv$obs - pred.cv$pred))# MAEsqrt(mean((pred.cv$obs - pred.cv$pred)^2))# RMSE## prediction ----predictionLocations <- readOGR(rpath, "predictionLocations")pred.discrete <- discretizePolygon(predictionLocations, cellsize = 1500, id = "ID")pred <- ataKriging(obs.discrete, pred.discrete, pointsv$pointVariogram) 2) 替代obs.discrete，并且拟合半方差函数 ## Replace the filen=18# 1-km数据av_link=paste('D://research//3_GHM//data//data_av_out_',n,'.csv',sep='')dp_link=paste('D://research//3_GHM//data//data_dp_out_',n,'.csv',sep='')# 原始数据：#av_link=paste('D://research//3_GHM//Result//2_av_dp//data_av_out_',n,'.csv',sep='')#dp_link=paste('D://research//3_GHM//Result//2_av_dp//data_dp_out_',n,'.csv',sep='')av=read.csv(av_link)dp=read.csv(dp_link)av=av[which(av$value>0),]dp=dp[which(dp$weight>0),]#av=read.csv('D://research//3_GHM//test//ata_1av.csv')#dp=read.csv('D://research//3_GHM//test//ata_1dp.csv')obs.discrete[["areaValues"]]=avobs.discrete[["discretePoints"]]=dpobs.discrete[["areaValues"]][["areaId"]]=as.character(obs.discrete[["areaValues"]][["areaId"]])obs.discrete[["discretePoints"]][["areaId"]]=as.character(obs.discrete[["discretePoints"]][["areaId"]])pointsv <- deconvPointVgm(obs.discrete, model="Exp", ngroup=10, rd=0.1, fig=TRUE)## cross validation ----pred.cv <- ataKriging.cv(obs.discrete, nfold=length(observations), pointsv)names(pred.cv)[6] <- "obs"pred.cv<-na.omit(pred.cv)summary(pred.cv[,c("obs","pred","var")])cor(pred.cv$obs, pred.cv$pred)# Pearson correlationmean(abs(pred.cv$obs - pred.cv$pred))# MAEsqrt(mean((pred.cv$obs - pred.cv$pred)^2))# RMSEpred.cv['re']=abs(pred.cv[4]-pred.cv[6])output_name=paste('pred_ghm_',n,'.csv',sep='')write.csv(pred.cv,output_name) 3）生成预测格网grid data ## prediction grid and zone attitudeimos=read.csv('D://research//3_GHM//Result//1_Spatial zones//spatial zones_2_1km.csv')coordinates(imos) <- c("lon","lat")#定义坐标 library(sf)library(sp)grd <- as.data.frame(spsample(imos, "regular", n=500))names(grd) <- c("X", "Y")coordinates(grd) <- c("X", "Y")# Add P's projection information to the empty gridproj4string(imos) <- proj4string(imos) # Temp fix until new proj env is adoptedproj4string(grd) <- proj4string(imos)## Spatial overlaya.data=over(grd,vector[,'zones'])grd$zone <- a.data$zonesgridded(grd) <- TRUE # Create SpatialPixel objectfullgrid(grd) <- TRUE # Create SpatialGrid object注意！grid数据的列名要修改成dp的格式。我太懒了，直接从Excel种修改。后续应该从代码中自己改好，这样全程都在R中进行。
这一步，千万要注意列名、x、y是否有问题。尤其是x和y！！！不要弄反。
predictionLocations <- readOGR(rpath, "predictionLocations")pred.discrete <- discretizePolygon(predictionLocations, cellsize = 15000, id = "ID")ap_grd=read.csv('D:/research/3_GHM/Result/3_prediction/grd.csv') ## 导入预测点数据pred.discrete$discretePoints=ap_grd # 把预测点数据赋值colnames(ap_grd) <- colnames(av) #av_grd_group=rbind(av,ap_grd) ##把预测点与原有area合并，生成新的area数据pred.discrete$areaValues=ap_grd # 赋值新area数据pred.discrete[["areaValues"]][["value"]]=NApred.discrete[["areaValues"]][["areaId"]]=as.character(pred.discrete[["areaValues"]][["areaId"]])pred.discrete[["discretePoints"]][["areaId"]]=as.character(pred.discrete[["discretePoints"]][["areaId"]]) 5) 预测格网的插值，并且可视化 #pred <- ataKriging(obs.discrete, pred.discrete, pointsv$pointVariogram,longlat = TRUE,showProgress=TRUE)pred <- ataKriging(obs.discrete, pred.discrete, pointsv$areaVariogram,longlat = TRUE,showProgress=TRUE)ggplot() + geom_sf() + coord_sf(ylim = c(151, 153),xlim=c(-34,-31), expand = FALSE) + theme_classic() + geom_point(data = pred, aes(x = centx, y = centy, color =pred ,alpha=0.5),scale=pred)+scale_colour_gradientn(colours = terrain.colors(10),breaks = c(-3,-2,-1,0,1,2,3),limits=c(-3,3))