Паншарпенинг при помощи 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) # получаем низкочастотную компоненту при помощи функции focal() пакета 'raster'
     
}

Обратите внимание, что мы оставили возможность задавать пользовательские значения для параметров filter и fun, что оставляет простор для возможных экспериментов.

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

processingPansharp <- function(pan, multi, filter = 'auto', fun = mean){
    ' Улучшает разрешение мультиспектрального растра за счёт панхроматического растра '
    # @param pan - растр с высоким разрешением, объект типа Raster
    # @param multi - растр с низким разрешением, чьё разрешение необходимо улучшить, объект типа Raster
    # @param filter - скользящее окно, объект типа матрица
    # @param fun - функция, которая производит вычисления на скользящем окне (параметр функции focal() )
    # @return pansharp - растер 'multi' с улучшенным разрешением, объект типа Raster
 
    LPF <- extractLPF(pan, multi, filter, fun)
         
    multi <- resample(multi, pan) # приводим экстент и разрешение мультиспектрального растра к панхроматическому
         
    all <- stack(multi, pan, LPF)
     
    bands <- nbands(multi)
    pan_band <- bands + 1
    lpf_band <- bands + 2
     
    # Производим паншарпенинг каналов мультиспектрального растра по очереди
    pansharp_bands <- list()
    for (band in 1:bands) {
        subset <- all[[c(band, pan_band, lpf_band)]]
        raster <- calc(subset, pansharpFun)
        pansharp_bands[[band]] <- raster
    }
     
    pansharp <- stack(pansharp_bands) # итоговый растр
}

Ну, и наконец, последняя функция, носящая чисто утилитарный характер (позволяет использовать минимум параметров при сохранении растров) - обёртка к функции сохранения растра из пакета 'raster'. Идея состоит в том, что мы почти всегда будем хотеть сохранить результат в формате GeoTIFF, и при этом иметь целочисленный тип данных. Дело в том, что в большинстве случаев значения пикселей мультиспектрального растра - целые числа, а после паншарпенинга они становятся числами с плавающей запятой. Десятичные значения пикселя не несут никакой полезной информации, а занимают лишние сотни мегабайт места на диске. Использование целых чисел в качестве типа данных растра экономит более 0,5 Гб дискового пространства при паншарпенинге сцен WorldView-2.

saveResult <- function(raster, path, format = 'GTiff', datatype = 'INT2S'){
    '  Сохраняет растр по указанному пути  '
    # @param raster - растр, который должен быть сохранён, объект типа Raser
    # @param path - путь к файлу (без расширения) - string
    # @param format - формат сохраняемого растра в соответствии с требованиями функции writeRaster() - string
    # @param datatype - тип данных сохраняемого растра в соответствии с требованиями функции to writeRaster() - string
     
    writeRaster(raster, 
                path, 
                format = format, 
                datatype = datatype, 
                overwrite = T)
}

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

Пахроматический растр с пространственным разрешением 0,5 м

Мультиспектральный растр с пространственным разрешением 2,0 м

Мультиспектральный растр после паншарпенинга в R

Результат можно сравнить с результатом паншарпенинга в Orfeo TollBox:

Мультиспектральный растр после паншарпенинга при помощи OTB

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