У меня ежедневные осадки с 61 гидропоста за 12 лет на водосборе (8000 км2).
Цель состоит в том, чтобы создать продукт суточного количества осадков с привязкой к сетке 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) для каждого дня отдельно. Однако для всех дней оба преобразования не подходят, как я тестировал с помощью теста Шапиро Вилка. Итак, я выбираю нормальное преобразование оценок (NCR), как было предложено Grimes & Pardo (2009) Геостатистический анализ осадков. и использовал проф. Код Эштон Шортридж
Приведенный ниже код используется для создания климатической вариограммы сезона дождей. Обратите внимание, я использовал дни, когда более 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 стандартной ошибкой (мм) или отклонением (мм2)?
Следует ли масштабировать климатологическую вариограмму (самородок и порог) с учетом дисперсии?
Обратное преобразование должно выполняться с отдельным параметром или без него?
Спасибо за помощь заранее.
