С координатами сильно промахнулись. Комментарии писал по ходу чтения, не правил (могли остаться лишние). В целом неплохо.
Код: Выделить всё
# R version 3.4.4 (2018-03-15) -- "Someone to Lean On"
library(geoR)
# --------------------------------------------------------------
# Analysis of Geostatistical Data
# For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
# geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
# --------------------------------------------------------------
# Руководствовался мануалом
# https://rpubs.com/Dr_Gurpreet/spatial_interpolation_kriging_R
fn = "T2M_DATA.csv"
tab = read.table(fn,header=TRUE,sep=",",as.is=TRUE)
tab2 <- tab[1:3]
str(tab2)
# 'data.frame': 35 obs. of 5 variables:
# $ LATI : int 559368 598947 638525 678104 717683 757263 796842 558854 598090 637325 ...
# $ LONG : int 4955453 4955939 4956669 4957642 4958858 4960319 4962023 5010996 5011483 5012212 ...
# $ T2M_aprel: num 11.99 11.03 8.8 5.59 3.02 ...
# преобразовываем в формат геоданных
tab2 <-as.geodata(tab2, coords.col = c(1, 2), data.col = 3)
str(tab2)
# List of 2
# $ LATI , LONG : int [1:35, 1:2] 559368 598947 638525 678104 717683 757263 796842 558854 598090 637325 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : NULL
# .. ..$ : chr [1:2] "LATI" "LONG"
# $ data : num [1:35] 11.99 11.03 8.8 5.59 3.02 ...
# - attr(*, "class")= chr "geodata"
summary(tab2)
# Оценка пространственной тенденции и нормальности набора данных
plot(tab2) # - не очень информативно
#---------------------------------------------------------------
data.lim = range(tab$T2M_aprel)
plot(LATI~LONG,data=tab,pch=20,cex=0.5+(tab$T2M_aprel-data.lim[1])/diff(data.lim)*2.5)
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# видим что-то похожее на тренд по диагонали, к тому же билинейный
# по хорошему, его нужно засовываить в универсальный кригинг, не знаю, есть ли он в этом пакете. будем смотреть
drift.model = lm(T2M_aprel~LATI*LONG,data=tab)
summary(drift.model)
#
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 6.055e+02 1.139e+02 5.314 8.72e-06 ***
# LATI -1.054e-03 1.677e-04 -6.286 5.45e-07 ***
# LONG -1.130e-04 2.248e-05 -5.027 1.99e-05 ***
# LATI:LONG 2.010e-10 3.309e-11 6.075 9.92e-07 ***
# Multiple R-squared: 0.907, Adjusted R-squared: 0.898
# видим, что тренд объясняет практически все (90% дисперсии).
# Вариограмму нужно считать для остатков от тренда
# Но пойдем по коду, будем искать, куда координаты пропали
tab$Resid = residuals(drift.model)
data.lim = range(tab$Resid)
plot(LATI~LONG,data=tab,pch=20,cex=0.5+(tab$Resid-data.lim[1])/diff(data.lim)*2.5)
hist(tab$Resid,col="gray")
qqnorm(tab$Resid,pch=20,cex=1.5)
qqline(tab$Resid,col="red")
ks.test(tab$Resid,"pnorm",mean(tab$Resid),sd(tab$Resid))
# остатки можно считать нормальными
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
plot(tab2, trend = "1st")
plot(variog(tab2, option = "cloud", max.dist = 161156.1))
# 161156.1 -- 1/2 максимальной дистанции
plot(variog(tab2, max.dist = 161156.1))
vario_model <- variog(tab2, trend = "1st", max.dist = 161156.1)
plot(vario_model)
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
vario_model.no_drift = variog(tab2, max.dist = 161156.1)
vario_model.drift = variog(tab2, max.dist = 161156.1, trend = "1st")
y.lim=c(0,max(vario_model.no_drift$v,vario_model.drift$v))
x.lim=c(0,161156.1)
plot(vario_model.no_drift$u,vario_model.no_drift$v,pch=20,col="black",cex=1.5,ylim=y.lim,xlim=x.lim)
points(vario_model.drift$u, vario_model.drift$v ,pch=20,col="red",cex=1.5)
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# Определение направленных эффектов/Анизотропии.
plot(variog4(tab2, trend = "1st", max.dist = 161156.1), omnidirectional = T)
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
vario_model.4dir = variog4(tab2, trend = "1st", max.dist = 161156.1)
# данных для определения анизотропии нет, там одни NA,
# поэтому рисуется только omnidirectional (последняя)
for(i in 1:5) print(as.integer(is.na(vario_model.4dir[[i]]$v)))
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# Подбор ковариационных моделей на вариограмме
exp_fit_fix <- variofit(vario_model, cov.model = "exponential", fix.nugget = T)
exp_fit_nofix <- variofit(vario_model, cov.model = "exponential", fix.nugget = F)
sph_fit_fix <- variofit(vario_model, cov.model = "spherical", fix.nugget = T)
sph_fit_nofix <- variofit(vario_model, cov.model = "spherical", fix.nugget = F)
plot(vario_model,pch=16)
lines(exp_fit_fix,col="green",lwd=4, lty = 1)
lines(exp_fit_nofix,col="red",lwd=4, lty = 2)
lines(sph_fit_fix,col="blue",lwd=4, lty = 3)
lines(sph_fit_nofix,col="black",lwd=4, lty = 4)
# Выбор наиболее подходящей модели.
# Полученная модель с наименьшим значением SSQ считается наиболее подходящей. Расчет SSQ для каждой подходящей модели:
(exp_SSQ_fix <- summary(exp_fit_fix)$sum.of.squares)
(exp_SSQ_nofix <- summary(exp_fit_nofix)$sum.of.squares)
(sph_SSQ_fix <- summary(sph_fit_fix)$sum.of.squares)
(sph_SSQ_nofix <- summary(sph_fit_nofix)$sum.of.squares)
which.min(list(exp_SSQ_fix,exp_SSQ_nofix, sph_SSQ_fix, sph_SSQ_nofix))
# оздание сетки прогнозов
prediction_grid <- expand.grid(seq(0, 800000, length.out = 100), seq(0, 800000, length.out = 100))
plot(prediction_grid)
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# смотрим, как наши точки лежат на сетке
points(LATI~LONG,data=tab,pch=20,cex=0.5+(tab$Resid-data.lim[1])/diff(data.lim)*2.5,col="red")
(tmp1=apply(prediction_grid,2,range))
# Var1 Var2
#[1,] 0e+00 0e+00
#[2,] 8e+05 8e+05
(tmp2=apply(tab[,1:2],2,range))
# LATI LONG
#[1,] 557285 4955453
#[2,] 796842 5184215
tmp2-tmp1
# LATI LONG
#[1,] 557285 4955453
#[2,] -3158 4384215
# сетка далеко в стороне от наших точек
prediction_grid <- expand.grid(seq(tmp2[1,2]-20000,tmp2[2,2]+20000,len=100),seq(tmp2[1,1]-20000,tmp2[2,1]+20000,len=100))
plot(prediction_grid)
points(LATI~LONG,data=tab,pch=20,cex=0.5+(tab$Resid-data.lim[1])/diff(data.lim)*2.5,col="red")
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# Кригинг
krig_temp <- krige.conv(tab2, loc=prediction_grid, krige=krige.control(obj.model=sph_fit_fix))
#gsph_fit_fix - смотрим по выводу команды which.min.
image(krig_temp,col=heat.colors(64))