如何反向转换正态分数转换后的数据

Question

我每天有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

作为日常雨量偏斜,变换i相测试立方根变换和日志变换(log1p)为每天单独.然而,对于所有这些日子来说,两种变换都不适合我用shapiro wilk测试测试.因此,我选择正常分数分数变换(NCR),如Grimes&Pardo(2009)降雨量地统计分析所示.并使用了Ashton Shortridge教授的代码

以下代码用于生成季风季节的气候变异函数.请注意,我曾使用超过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

类似的kriging读取每天非零,我已经规模每天下雨并改变它.不确定以下方法是否正确？

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)

我的困惑如下:

1. var1.var是标准误差(mm)还是方差(mm2)？

2. 气候变异函数(金块和基石)是否应该根据方差进行缩放？

3. 后转换应该单独使用还是不选择？

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

Answer 1

AE Journel（数学地质学，大约 1980 年）详细讨论了对数正态变换数据的克里金法与对数正态克里金法。它显示了在克里金法之前使用数据的任何非线性变换的困难。请注意，该期刊已更名多次，现在是Mathematical Geosciences

对于降雨量，您还需要考虑数据是否代表一段时间内（一小时、一天、一周等）的累积量

如果感兴趣区域的海拔变化很大，那么您需要合并数据位置的海拔。在文献中至少有两种方法可以做到这一点；(1) 将降雨平均值建模为函数