我每天有61个测量站的降雨量，为期12年（8000平方公里） .

目标是创造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

由于每日降雨量有所偏差，我将使用 cube root transform 和 log transform (log1p) 进行转换测试 each day individually . 但是对于所有这些日子来说，这两个变换都不适合我用shapiro wilk测试测试 . 因此，我选择 Grimes & Pardo(2009) Geostatistical analysis of rainfall 建议的正常分数分数变换（NCR） . 并使用了Prof. Ashton Shortridge code

以下代码用于生成季风季节的气候变异函数 . 请注意，我曾使用超过30％的电台报告下雨的日子 . 没有找到任何参考 . 当我达到65％的天数超过阈值时，选择30％ .

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

Similary for kriging read each day non-zero, i had scaleed daily rain and transform it . 不确定以下方法是否正确？

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）？

  • 气候变异函数（金块和窗台）是否应根据方差进行缩放？

  • 后退转换应该单独使用还是不选择？

我在这里先向您的帮助表示感谢 .

r interpolation kriging gstat