Помогите разобраться как построить именно карту различий между растрами с использованием Time-Weighted Dynamic Time Warpin?
Сами различия между растрами, образующие временной ряд определил с помощю библиотек "twdtw" и dtw.
Код: Выделить всё
library(terra)
library("twdtw")
library(dtw)
# Установка рабочей директории
Data <- "/home/taras/GRIDS"
setwd(Data)
list.files(pattern="\\.tif$")
# Список файлов
raster_files <- list.files(pattern = "^summer_\\d{4}\\.tif$", full.names = TRUE)
# Загрузка растров
time_series <- rast(raster_files)
# Извлечение годов из имен файлов
years <- as.numeric(gsub("summer_(\\d{4})\\.tif", "\\1", basename(raster_files)))
# Присвоение временных меток
time(time_series) <- as.Date(paste0(years, "-06-15")) # 15 июня как пример
# Проверка результата
print(time(time_series))
# Эталонная временная серия (в данном случае синусоида как пример)
ref_times <- seq(as.Date("1984-06-15"), as.Date("2023-06-15"), by = "year")
ref_values <- sin(2 * pi * seq_along(ref_times) / 6) # Синусоида
reference <- data.frame(date = ref_times, value = ref_values)
# Пример динамического выравнивания для каждого растрового временного ряда
distances <- numeric(length(raster_files))
for (i in 1:length(raster_files)) {
# Извлечение данных для каждого растрового временного ряда
raster_values <- values(time_series[[i]])
# Преобразование векторных значений растров в одномерный вектор
raster_values <- as.vector(raster_values)
# Убедитесь, что длина временных рядов одинаковая, если необходимо, интерполируйте
if(length(raster_values) != length(ref_values)) {
# Пример: интерполяция значений растров на длину эталона
raster_values <- approx(raster_values, n = length(ref_values))$y
}
# Вычисление DTW расстояния
alignment <- dtw(raster_values, ref_values, keep = TRUE)
# Сохранение расстояния DTW для каждого растрового временного ряда
distances[i] <- alignment$distance
}
# Вывод результатов
data.frame(year = years, dtw_distance = distances)
Линейный график DTW расстояния по годам:
library(ggplot2)
# Создание data frame
results <- data.frame(year = years, dtw_distance = distances)
# Построение линейного графика
ggplot(results, aes(x = year, y = dtw_distance)) +
geom_line(color = "blue") +
geom_point(color = "red") +
labs(title = "DTW Distance Over Years",
x = "Year",
y = "DTW Distance") +
theme_minimal()
Гистограмма распределения DTW расстояний:
# Гистограмма распределения DTW расстояний
ggplot(results, aes(x = dtw_distance)) +
geom_histogram(binwidth = 0.1, fill = "lightblue", color = "black") +
labs(title = "Distribution of DTW Distances",
x = "DTW Distance",
y = "Frequency") +
theme_minimal()Попытка в лоб найти различия приводят к результату на рисунке
Код: Выделить всё
# Загрузка растров для 1984 и 2023 года
raster_1984 <- rast("summer_1984.tif")
raster_2023 <- rast("summer_2023.tif")
# Определение размеров блока (например, 100x100 пикселей)
block_size <- 100
# Получаем размеры растров
ncols <- ncol(raster_1984)
nrows <- nrow(raster_1984)
# Создание пустого объекта для хранения результата
difference_raster <- raster_1984
values(difference_raster) <- NA # Инициализация пустыми значениями
# Обрабатываем блоки по частям
for (i in seq(1, nrows, block_size)) {
for (j in seq(1, ncols, block_size)) {
# Ограничиваем индексы на случай, если блок выходит за пределы
row_idx <- i:min(i + block_size - 1, nrows)
col_idx <- j:min(j + block_size - 1, ncols)
# Создаем координаты для текущего блока
grid <- expand.grid(row = row_idx, col = col_idx)
# Извлекаем значения для текущего блока
block_1984 <- extract(raster_1984, grid)
block_2023 <- extract(raster_2023, grid)
# Диагностика: выводим значения для блока
print(paste("Block 1984:", toString(block_1984)))
print(paste("Block 2023:", toString(block_2023)))
# Убираем NA значения
valid_indices <- !is.na(block_1984) & !is.na(block_2023)
block_1984 <- block_1984[valid_indices]
block_2023 <- block_2023[valid_indices]
# Диагностика: выводим очищенные блоки
print(paste("Valid Block 1984:", toString(block_1984)))
print(paste("Valid Block 2023:", toString(block_2023)))
# Применяем DTW для текущего блока
if (length(block_1984) > 0 && length(block_2023) > 0) {
dtw_result <- dtw(block_1984, block_2023)
# Диагностика: выводим результат DTW
print(paste("DTW Distance:", dtw_result$distance))
# Сохраняем результат в соответствующем месте
diff_results[row_idx, col_idx] <- dtw_result$distance
}
}
}
# Преобразуем результаты в новый растр
raster_diff <- raster_1984
values(raster_diff) <- diff_results
# Сохранение нового растра
writeRaster(raster_diff, "difference_raster.tif", overwrite=TRUE)
# Визуализация разницы
plot(difference_raster, main = "DTW Difference between 2023 and 1984 (Block-wise)", col = terrain.colors(100)) 
