Паншарпенинг при помощи R

Эта страница является черновиком статьи.

В данной статье рассказывается о том как произвести паншарпенинг в программной среде R. [via Misanthrope's Thoughts]

Введение

Одной из наиболее рапространённых задач при работе со спутниковыми снимками является паншарпенинг. С этой задачёй успешно справляются и проприетарные и открытые программы для работы с ДДЗЗ. Но в некоторых случаях вам может захотеться иметь больший контроль над проводиммыми операциями или, быть может использовать алгоритм, ещё не реализованный ни в одном ПО. В таком случае на помощь придёт один из языков программирования, например, R.

В данной статье будет продемонстрирована реализация наиболее простого алгоритма паншарпенига - модуляции высоких частот.

Немного теории

На удивление сложно найти определение для слова "паншарпенинг". Попробуем сформулировать его своими словами, взяв за основу описание из книги Шовергердт Р.А. Дистанционное зондирование. Модели и методы обработки изображений - М.: Техносфера, 2010. (с. 416):

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

Одним из наиболее простых алгоритмов паншарпенинга является модуляция высоких частот. В этом алгоритме для получения улучшенного изображения в канале PXS панхроматический снимок PAN попиксельно умножается на изображение с низким разрешением XS, после чего полученный результат нормируется на низкочастотную компоненту панхроматического снимка PAN_smooth (панхроматический снимок, обработанный сглаживающим фильтром с окном, соответствующим размеру пикселя изображения с низким разрешением):

где i и j - индексы пикселей.

Реализация

Для работы с растрами, имеющими географическую привязку в R есть замечательный пакет raster. На основании функций, доступных в этом пакете мы и разработаем наш скрипт. В использовании он будет чрезвычайно прост:

library(raster)

pan <- raster('pan.tif') # загружаем панхроматический растр
multi <- brick('multi.tif') # загружаем мультиспектральный растр
pansharp <- processingPansharp(pan, multi) # производим паншарпенинг (эту функцию мы создадим)
output_path <- 'путь_и_имя_файла_без_расширения' # зададим путь для сохранения результата
saveResult(pansharp, output_path) # сохраним результат

Начнём создавать необходимые функции. Первая функция - самая простая. Она всего-лишь реализовывает формулу, представленную выше.

pansharpFun <- function(raster){
    ' Эта функция производит паншарпенинг  ' 
    # @param raster - объект типа Raster с тремя каналами: 1)с низким разрешением, 2) с высоким разрешнеием, 3) низкочастотная компонента растра с высоким разрешением
    # @return pansharpened_raster - подвергшийся паншарпенингу  объект типа Raster
    # pansharp = Lowres * Highres / LPF[Highres]
     
    pansharpened_raster <- (raster[,1] * raster[,2]) / raster[,3]
}

Вторая функция будет вычислять низкочастотную сотавляющую растра с высоким расзрешением.

extractLPF <- function(pan, multi, filter = 'auto', fun = mean) {
    ' Возвращает низкочастотную компоненту растра с высоким разрешением 
        с использованием скользящего окна, соответствующего размеру пикселя растра с низким разрешением '
    # @param pan - растр с высоким разрешением, объект типа Raster
    # @param multi - растр с низким разрешением, который в последствии будет подвергнут паншарпенингу, объект типа Raster
    # @param filter - скользящее окно, объект типа матрица
    # @param fun - функция, которая производит вычисления на скользящем окне (параметр функции focal() )
    # @return LPF - низкочатотная компонента растра с высоким разрешением, объект типа Raster
     
    # Расчитаем размер скользящего окна, если его параметры не заданы
    if (filter == 'auto') {
        pan_res <- res(pan) # (x, y) пространственное разрешение растра с высоким разрешением в единицах CRS
        multi_res <- res(multi) # (x, y) пространственное разрешение растра с низким разрешением в единицах CRS
        x_res_ratio <- round(multi_res[1]/pan_res[1])
        y_res_ratio <- round(multi_res[2]/pan_res[2])
        total <- x_res_ratio + y_res_ratio
        filter <- matrix(1, nc = x_res_ratio, nr = y_res_ratio)
        }
         
        # Убедимся, что у скользящего окна нечётное количество столбцов и строк (требование функции focal() )
        if (nrow(filter)%%2 == 0) {
            filter <- rbind(filter, 0)
        } 
        if (ncol(filter)%%2 == 0) {
            filter <- cbind(filter, 0)
       
         
    LPF <- focal(pan, w = filter, fun = fun) # получаем низкочастотную компоненту
     
}

Результаты и обсуждение

К сожалению, в raster не реализована многопоточность, поэтому производительность операций по обработке больших растров оставляет желать лучшего.