Кригинг, теряются координаты.

[_taras_]

Доброго времени!
Пробую разобраться с кригингом и возникли вопросы.
1 - где-то в процессе или теряются или неверно читаются координаты точек.
2 - как сохранить полученный результат в необходимой системе координат.
Данные представляют таблицу с колонками -
LAT, LON, Jan
44.75,27.75,1.0
Координаты в десятичных градусах, epsg 4326
Преобразовываю в формат геоданных

Код: Выделить всё

tab2 <-as.geodata(tab2, coords.col = c(1, 2), data.col = 3)

Выбираю подходящую модель, создаю сетку прогнозов, делаю пространственное прогнозирование и смотрю результат. По мануалу результат привязан к координатам. У меня на рисунке

kriging.jpeg (63.38 КБ) 3980 просмотров

В реальности должны быть на точках

points.jpeg (189.32 КБ) 3980 просмотров

[gamm]

1. Какая программа/пакет (судя по всему, в среде R). Приведите весь скрипт и данные, можно модельные. Там кригингов много разных.

2. Интерполировать в градусах очень плохая идея. Надо сначала точки перепроецировать

[_taras_]

Выкладываю
T2M_DATA.csv - сами данные. Колонки LATI, LONG - координаты для WGS_1984_UTM_ZONE_35N [+proj=utm +a=6378137.0 +b=6356752.314245 +zone=35 +no_def]
LAT, LON - WGS_1984, epsg 4326
Kriging_R.txt - что делаю.

[gamm]

Выкладываю

С координатами сильно промахнулись. Комментарии писал по ходу чтения, не правил (могли остаться лишние). В целом неплохо.

Код: Выделить всё

# 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))

[gamm]

строго говоря, приведенный выше метод определения параметров вариограммы имеет ряд больших недостатков:

1) проблема "курица или яйцо" - для вычисления вариограммы нужны остатки от тренда, а для правильного вычисления тренда нужна вариаограмма. Круг замкнулся.

2) способ построения модели вариограммы основан на вычислении среднего в секторах (по дальности и направлению). С этим проблема - квадраты разности дают выбросы, размер секторов и прочее непонятно как выбирать, и разный выбор дает сильно разную картину. Точность подгонки модели к вычисленным значениям роли не играет, поскольку сами значения имеют мало общего с реальностью.

Наука статистика им.тов.Cressie предполагает совместное определение тренда и параметров вариограммы без всякого вычисления эксперименьальной вариограммы, за счет построения общей функции правдоподобия и ее максимизации. Численный вариант такого подхода есть в R, это gls() в пакете nlme, и (для нелинейного тренда) gamm() в пакете mgcv. Ниже пример определения параметров вариограммы вместе с трендом, и выбор наилучшей модели.

Там еще много чего можно, но об этом в другой раз

Код: Выделить всё

#=========================================================================
# --- variogram model with drift
#=========================================================================
library(MASS) 
library(nlme)

lm3 <-gls(T2M_aprel~LATI+LONG,data=tab,method="ML")
lm3a <-gls(T2M_aprel~LATI*LONG,data=tab,method="ML")
anova(lm3,lm3a)

lm3 <-gls(T2M_aprel~LATI*LONG,data=tab,method="REML")

# модели с различными вариантами корреляции, метод REML
lm3c<-gls(T2M_aprel~LATI*LONG,data=tab,method="REML",correlation = corLin(nugget=TRUE,form = ~ LATI+LONG))
lm3d<-gls(T2M_aprel~LATI*LONG,data=tab,method="REML",correlation = corSpher(nugget=TRUE,form = ~ LATI+LONG))

# выбираем корреляцию
anova(lm3,lm3c)
anova(lm3,lm3d)

coef(lm3d$modelStruct$corStruct, unconstrained = FALSE)
#        range       nugget 
# 1.696779e+05 3.697991e-10 

[_taras_]

Код: Выделить всё

gamm, спасибо огромное за Ваши комментарии.
С координатами оказалось всё просто... Перепутаны широта и долгота. Исправил и точки на месте.
Остался вопрос как сохранить результат кригинга в geoTIF?

Код: Выделить всё

krig_temp <- krige.conv(tab2, loc=prediction_grid, krige=krige.control(obj.model=sph_fit_fix))
image(krig_temp,col=heat.colors(64))

У меня задача составить на основании имеющихся данных по температурам (воздуха и поверхности) их градиент и посмотреть насколько она отличается для нескольких природно-заповедных территорий, находящихся в этом регионе.

[gamm]

Код: Выделить всё

library(raster)

(tmp2=apply(tab[,1:2],2,range))

# многие программу требуют квадратных пикселей
step=min(diff(range(tmp2[,1]))/100,diff(range(tmp2[,2]))/100)

x.coord = seq(tmp2[1,1]-20000,tmp2[2,1]+20000,by=step)
y.coord = seq(tmp2[1,2]-20000,tmp2[2,2]+20000,by=step)
nx=length(x.coord); xmin=min(x.coord); xmax=max(x.coord)
ny=length(y.coord); ymin=min(y.coord); ymax=max(y.coord)

prediction_grid <- expand.grid(x.coord,y.coord)
plot(prediction_grid)
points(LONG~LATI,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))

krige_raster <- raster(nrows=ny,ncols=nx,xmn=xmin,xmx=xmax,ymn=ymin,ymx=ymax,crs=prj.UTM35,vals=krig_temp$predict)
krige_raster[,] = krige_raster[nrow(krige_raster):1,]

plot(krige_raster,col=gray(0:64/64))
image(krig_temp,col=gray(0:64/64))

# Сохраните растр в формате GeoTIFF
writeRaster(krige_raster, filename="krige_result.tif", format="GTiff", overwrite=TRUE)