Задать очграничение в коде для рассчетов

[_taras_]

Доброго времени!
Нужно посчитать индекс для списка растров в которых есть пиксели = 0 (бывшие облака). Если сам подсчет в целом понятен, то как прописать, то что пиксели с 0 значениями считать не надо и присвоить NA?
Мой код

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

library(terra)
library(sf)
rm(list = ls())

# Устанавливаем рабочую директорию
Data <- "/home/taras/GRIDS"
setwd(Data)

# Загружаем векторный слой с границами
dnestr_area <- vect("~/GRIDS/Dnestr_area/Dnestr_area.shp")

# Создаем папку для выходных файлов, если её нет
output_dir <- file.path(Data, "TCI_results")
if (!dir.exists(output_dir)) {
  dir.create(output_dir)
}

# Получаем список растровых файлов
rast_list <- list.files(pattern = "TEMP_\\d{4}_\\d{2}_15\\.tiff$", full.names = TRUE)

# Обрабатываем каждый файл
for (file in rast_list) {
  # Загружаем растр
  TEMP_rast <- rast(file)
  
  # Проверяем, совпадает ли проекция с векторным слоем, если нет - преобразуем
  if (!same.crs(TEMP_rast, dnestr_area)) {
    dnestr_area <- project(dnestr_area, crs(TEMP_rast))
  }
  
  # Обрезаем по границам полигона
  TEMP_crop <- crop(TEMP_rast, dnestr_area)
  TEMP_mask <- mask(TEMP_crop, dnestr_area)
  
  # Вычисляем TEMPmin и TEMPmax в пределах полигона
  TEMP_min <- -53.6567573547363
  TEMP_max <- 85.494255065918
  
  # Проверяем деление на ноль
  if (is.na(TEMP_min) || is.na(TEMP_max) || (TEMP_max - TEMP_min) == 0) {
    warning(paste("Пропущен файл", file, "- некорректные TEMPmin/TEMPmax"))
    next
  }
  
  # Вычисляем TCI
 TCI <- (TEMP_max - TEMP_mask)/ (TEMP_max - TEMP_min)

  # Формируем имя выходного файла
  TCI_filename <- sub("TEMP_", "TCI_", basename(file))
  TCI_path <- file.path(output_dir, TCI_filename)
  
  # Сохраняем результат
  writeRaster(TCI, filename = TCI_path, overwrite = TRUE)
  
  cat("Обработан:", TCI_filename, "\n")
}

cat("Готово! Все файлы сохранены в:", output_dir, "\n")

[Константин Силкин]

Здравствуйте! Вы там что-то пропустили? Сперва присваиваете константы переменным TEMP_min и TEMP_max, а затем проверяете их значение. Не понял смысла

[_taras_]

Здравствуйте! Вы там что-то пропустили? Сперва присваиваете константы переменным TEMP_min и TEMP_max, а затем проверяете их значение. Не понял смысла

Спасибо за участие.
Это неудаленный кусок, проверяющий корректность конкретных TEMPmin и TEMPmax полученных после

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

  TEMP_min <- as.numeric(global(TEMP_mask, fun = "min", na.rm = TRUE))
  TEMP_max <- as.numeric(global(TEMP_mask, fun = "max", na.rm = TRUE))  

Это ещё один вопрос требующий элегантного решения -- как найти экстремальные значения TEMPmin и TEMPmax в списке растров. Т.е. имеется ..надцать растров и надобно найти самое минимальное (-15,25) и самое максимальное (45,53) значение. Сейчас решил в лоб отдельным скриптом.

[Константин Силкин]

Тогда что в итоге-то надо решить?

[_taras_]

Тогда что в итоге-то надо решить?

Как сделать так, чтобы пиксели с определенными значениями не учитывались? Т.е. значение пикселя = 0 формула не применяется и в результриующий растр для этого пикселя пишется 0 или NA

[Константин Силкин]

Сделайте индекс. Не в смысле вашего ДЗЗ-индекса, а в смысле индексации элементов массива. Применяйте этот индекс как в левой, так и в правой части оператора присваивания

[_taras_]

Сделайте индекс. Не в смысле вашего ДЗЗ-индекса, а в смысле индексации элементов массива. Применяйте этот индекс как в левой, так и в правой части оператора присваивания

Не очень понял, что Вы предлагаете, но натолкнули на идею использовать маски туч (как чувствовал, что пригодиться)... Просто перемножу и готовый результат...
Но хотел бы увидеть Вашу реализацию в качестве примера, поскольку учусь в процессе...

[Константин Силкин]

По большому счёту маска в этом смысле – это и есть тот индекс, что я предлагаю. Сделайте сами маску нулевых значений, маскируйте входные и выходные массивы и производите спокойно свои вычисления.
Простите, на R не пишу уже 4 года, поэтому могу спутать мелкие нюансы синтаксиса, так что писать код не стану.
Предлагаю только логику, которую я сам всегда применяю, на каком языке это бы не писалось

[gamm]

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

library(raster)
foo=raster("foo.tif")
foo[foo[]==0]=NA

[_taras_]

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

library(raster)
foo=raster("foo.tif")
foo[foo[]==0]=NA

В своем коде добавил после

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

  # Вычисляем TCI
 TCI <- (TEMP_max - TEMP_mask)/ (TEMP_max - TEMP_min)
 # Замена значений TCI на 0 для пикселей, где TEMP_mask = 0
 TCI[TEMP_mask == 0] <- NA

И вроде все отлично...

[gamm]

я не очень понимаю поверх чего написан пакет terra, но для пакета raster есть принципиальная разница между foo==0 и foo[]==0, особенно для больших файлов (по скорости и памяти): в первом случае все делается на диске, через небольшой буфер, во втором - растр целиком грузится в память, что может резко ускорить работу (если памяти хватает). Поэтому я предпочитаю грузить все в память, и в ней считать, оперируя с векторами и матрицами, а потом положить все на место и сохранить. С полигонами аналогично, он переводится в номера ячеек (аналог маски), и с этим подмножеством можно работать.