我在一个流域（8000 平方公里）的 12 年内每天都有来自 61 个测量站的降雨量。
目标是创建 5 公里和 25 公里分辨率的网格化每日降雨量产品。由于站点数量很少，即使在雨季也不是所有站点都下雨，所以我选择使用气候变异函数。
典型一天（irain）的雨量计测量值如下，NA 表示的缺失值很少。
7.8 4.4 15.4 19.1 5.8 0 42 6.4 21 21 0 0 0 15.6 0 0 10 5 1.2 0 14.4 NA 25 13.2 0 9.2 2 6.6 7.8 13.2 15.4 NA 9 0 15.5 0 18.6 6 0 4.8 10.6 0 9 0 12.4 NA 12 0 3 14 10 10 0 68 21.8 18 14.8 5.4 7 0 NA
由于每日降雨量存在偏差，为了转换，我分别用立方根变换和对数变换 (log1p)对每一天进行了测试。然而，在我用 shapiro wilk 测试进行测试时，这两种转换都不适合。因此，我按照Grimes & Pardo(2009) Geostatistical analysis of rains的建议选择了正态分数转换（NCR）。并使用了 Ashton Shortridge 教授的代码
下面的代码用于生成季风季节的气候变异函数。请注意，我使用了超过 30% 的站点报告下雨的日子。没有找到任何参考资料。选择 30%，因为我有大约 65% 的天数可以通过阈值。
lag = 3
bins.vario = seq(0, 75, lag)
nb.bins = length(bins.vario) - 1
nb.classes = numeric(nb.bins)
vario.emp = array(0,c(nb.bins,6))
variance.emp = array(NA,c(10000,nb.bins))
vario.emp = as.data.frame(vario.emp)
class(vario.emp) = c("gstatVariogram","data.frame")
names(vario.emp) = c("np","dist","gamma","dir.hor","dir.ver","id")
nRows = nrow(stn.rain.subset)
for (i in 1:nRows)
{
irain = stn.rain.subset[i,]
isMissing = is.na(irain)
isZero = (irain == 0)
irain = irain[!isMissing & !isZero]
irain = as.numeric(irain)
rain.mean[[i]] = mean(irain); rain.var[[i]] = var(irain);
# Testing with log-transformation
# irain.logtt = log1p(irain)
# # qqPlot(irain.logtt,distribution = "norm")
# res = shapiro.test(irain.logtt)
# pval[[i]] = res$p.value
if (var(irain)>0)
{
# Scaling of the rain on each day. rain in mm. After scaling this is unitless
irain.scaled = irain/sqrt(var(irain))
irain.nsc = nscore(irain)
score = irain.nsc$nscore
easting = lon.UTM[!isMissing & !isZero] # Removing the stn.locs with NA values
northing = lat.UTM[!isMissing & !isZero]
rain = data.frame(easting,northing,score)
names(rain) = c("easting","northing","rain")
coordinates(rain) = c("easting","northing")
proj4string(rain) = UTM
v = variogram(rain~1, data = rain,boundaries = bins.vario)
bin.nb = ceiling(v$dist/lag)
nb.classes[bin.nb] = nb.classes[bin.nb]+1
vario.emp[bin.nb,] = vario.emp[bin.nb,]+v
}
非零降雨的气候变异函数：
同样，我生成了指标变异函数。
现在的问题是如何使用气候变异函数进行克里金法。
"model" "psill" "range" "kappa" "ang1" "ang2" "ang3" "anis1" "anis2"
"Nug" 0.446609415762287 0 0 0 0 0 1 1
"Sph" 0.533499909345213 51.7206027419321 0.5 0 0 0 1 1
与克里金法类似，每天读取非零，我已经缩放了每天的降雨量并对其进行了转换。不确定以下方法是否正确？
rain = data.frame(easting,northing,score)
# Multiplying the nugget and sill with var(rain) for each day.
clim.vrmod$psill = clim.vrmod$psill * var(irain)
krige.ok = krige(rain[,3]~1,locations =~easting+northing,data = rain,newdata=output.grd,model = clim.vrmod)
krige.ok$var1.pred.bt = (backtr(krige.ok$var1.pred,irain.nsc, tails='separate'))*sqrt(var(irain))
krige.ok$var1.se = (krige.ok$var1.var)
我的困惑如下：
var1.var 是标准误差（mm）还是方差（mm2）？
气候变异函数（块金和窗台）是否应该随方差缩放？
应该使用单独或无选项进行反向转换？
我在这里先向您的帮助表示感谢。