library(sp)
library(geoR)
library(rgeos)
library(raster)

# Чтение данных
tab <- read.csv(file="/home/taras/test.csv", header=TRUE, sep=",", dec=".")
# Выбираем нужные колонки (например, c 3-й по n-ю колонку для обработки)
columns_to_process <- colnames(tab)[3:ncol(tab)]  # Замените '3' на индекс первой колонки с данными для обработки
# Цикл по колонкам
for (col in columns_to_process) {
  
  # Преобразуем данные в формат geodata, используя текущую колонку
  tab2 <- as.geodata(tab, coords.col = c(1, 2), data.col = which(colnames(tab) == col))
  
  # Модель дрейфа
  drift.model = lm(paste0(col, " ~ LAT*LONG"), data=tab)
  tab$Resid = residuals(drift.model)
   # Модель дрейфа
  drift.model = lm(paste0(col, " ~ LAT*LONG"), data=tab)
  tab$Resid = residuals(drift.model)
  data.lim = range(tab$Resid)
  
  # Вычисление вариограммы
  vario_model <- variog(tab2, trend = "1st", max.dist = 161156.1)
  
  # Подгонка модели сферической вариограммы
  sph_fit_fix <- variofit(vario_model, cov.model = "spherical", fix.nugget = TRUE)

  # модель
(sph_SSQ_fix <- summary(sph_fit_fix)$sum.of.squares)
  
  
  # Создание сетки прогнозов
  prediction_grid <- expand.grid(seq(4844137, 5299450, length.out = 900), 
                                                    seq(479872, 882442, length.out = 900))
  
  # Кригинг
  krig_temp <- krige.conv(tab2, loc=prediction_grid, krige=krige.control(obj.model=sph_fit_fix))
  krige <- krige.control(obj.model=sph_fit_fix))
 
  (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)

  krig_temp <- krige.conv(tab2, 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)
  
  # Переворот растра для сохранения
  krige_raster[,] = krige_raster[nrow(krige_raster):1,]
  
  # Генерация имени файла для сохранения
  output_filename <- paste0("/home/taras/GRIDS/kriging_result_", col, ".tif")
  
  # Сохранение растра в GeoTIFF
  writeRaster(krige_raster, filename=output_filename, format="GTiff", overwrite=TRUE)
  
  print(paste("Сохранен результат для колонки:", col))
}