【发布时间】：2017-10-06 14:29:33
【问题描述】：
我在一个集水区（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 rain 建议的正态分数转换（NCR）。并使用Prof. Ashton Shortridge code
以下代码用于生成季风季节的气候变异函数。请注意，我使用了超过 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）？
是否应缩放气候变异函数（块金和窗台） 有方差吗？
应该使用单独或无选项进行反向转换？
提前感谢您的帮助。
【问题讨论】：
-
你需要把你的问题分解成碎片。否则你不会得到答案。我在同一领域工作，经过几行之后，您的问题对我来说变得势不可挡，因为您的问题是针对您的，并且更像是您在研究中面临的问题。期望每个人都通过几行文字明白这一点是不现实的。
-
@马苏德。感谢您提出我提出的问题。作为一名研究人员，我们觉得我尽了最大努力，并提出了一个尽可能简单的问题，以及为什么它看起来很复杂。我会将其分解为 2 个或更多问题。
-
实际上，如果您给我们一个示例，说明您想做什么以及您尝试了什么，而不必依赖背后的理论，那将会有所帮助。但是您的问题似乎集中在理论上，因此您最好在 stats.stackexchange.com 上尝试
标签： r interpolation kriging gstat