2) исправил методические ошибки с трендом на скорую руку.
3) не хватает анализа остатков, как минимум на нормальность, гетероскедастичность и пр.
Код: Выделить всё
library(sp)
library(geoR)
library(rgeos)
library(raster)
library(MASS)
library(mgcv)
dat.dir = "d:/temp/test"
fn=sprintf("%s/test.csv",dat.dir)
p.krig = read.table(fn,header=TRUE,sep=";")
p.krig = p.krig[,c(2:1,3:ncol(p.krig))] # иначе картинка вывернута, LONG->x, LAT->y
plot(p.krig[,1:2],pch=20,cex=(p.krig[,3]-min(p.krig[,3]))/diff(range(p.krig[,3]))*2)
hist(p.krig[,3])
# [1] "LAT" "LONG" "y_84_01" "y_84_02" "y_84_03"
hist(dist(p.krig[,1:2]))
M0a = lm(y_84_01~LONG+LAT, data=p.krig)
M0b = lm(y_84_01~LONG*LAT, data=p.krig)
anova(M0a,M0b)
# тренд билинейный
M0 = gls(y_84_01~LONG*LAT, data=p.krig)
M1 = gls(y_84_01~LONG*LAT, data=p.krig, correlation=corLin (form = ~ LONG + LAT,nugget=TRUE))
M2 = gls(y_84_01~LONG*LAT, data=p.krig, correlation=corSpher(form = ~ LONG + LAT,nugget=TRUE))
M3 = gls(y_84_01~LONG*LAT, data=p.krig, correlation=corExp (form = ~ LONG + LAT,nugget=TRUE))
anova(M0,M1,M2,M3)
# сферическая фариограмма лучшая, nuget не нужен.
summary(M2)
# Parameter estimate(s):
# range nugget
#2.423953e+05 1.529560e-10
# максимальное расстояние в районе 3e5 - 4e5
max_dist = 2.5e5
# Чтение данных
tab <- p.krig
# Выбираем нужные колонки (например, c 3-й по n-ю колонку для обработки)
columns_to_process <- colnames(tab)[3:ncol(tab)] # Замените '3' на индекс первой колонки с данными для обработки
# Цикл по колонкам
col = columns_to_process[1]
# Преобразуем данные в формат geodata, используя текущую колонку
tab2 <- as.geodata(tab, coords.col = c(1, 2), data.col = which(colnames(tab) == col))
tab$Resid = residuals(M2,type="response")
# Вычисление вариограммы
vario_model <- variog(tab2, data=tab$Resid, max.dist = max_dist)
plot(vario_model)
# Подгонка модели сферической вариограммы
sph_fit_fix <- variofit(vario_model, cov.model = "spherical", fix.nugget = TRUE, nugget = 0)
# модель
(sph_SSQ_fix <- summary(sph_fit_fix)$sum.of.squares)
(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],tmp2[2,1],by=step)
y.coord = seq(tmp2[1,2],tmp2[2,2],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(p.krig[,1:2],pch=20,cex=(p.krig[,3]-min(p.krig[,3]))/diff(range(p.krig[,3]))*2,col="red")
# чисто посмотреть, в каком порядке укладываются значения
krig_temp <- krige.conv(tab2, loc=prediction_grid, krige=krige.control(obj.model=sph_fit_fix))
r_test <- raster(nrows=ny, ncols=nx, xmn=xmin, xmx=xmax, ymn=ymin, ymx=ymax, crs=crs.geo, vals=krig_temp$predict)
plot(r_test)
points(p.krig[,1:2],pch=20,cex=(p.krig[,3]-min(p.krig[,3]))/diff(range(p.krig[,3]))*2)
# кригинг остатков - как засунуть в него билинейный тренд не нашел
krig_temp <- krige.conv(tab2, data=tab$Resid, loc=prediction_grid, krige=krige.control(obj.model=sph_fit_fix))
crs.geo <- CRS("+proj=utm +zone=35 +datum=WGS84 +units=m +no_defs")
krige_raster <- raster(nrows=ny, ncols=nx, xmn=xmin, xmx=xmax, ymn=ymin, ymx=ymax, crs=crs.geo, vals=krig_temp$predict)
plot(krige_raster)
# модель тренда
trend.val = predict(M2,newdata=data.frame(LONG=prediction_grid[,1],LAT=prediction_grid[,2]),type="response")
trend_raster = krige_raster
trend_raster[] = 0
trend_raster[,] = matrix(trend.val,ncol=nx,nrow=ny,byrow=TRUE)
plot(trend_raster)
# сложить тренд и модель остатков
krige_raster[] = krige_raster[] + trend_raster[]
plot(krige_raster)
# Переворот растра для сохранения
krige_raster[,] = krige_raster[nrow(krige_raster):1,]
plot(krige_raster)
# Генерация имени файла для сохранения
output_filename <- paste0("/home/taras/GRIDS/kriging_result_", col, ".tif")
# Сохранение растра в GeoTIFF
writeRaster(krige_raster, filename=output_filename, format="GTiff", overwrite=TRUE)
print(paste("Сохранен результат для колонки:", col))
