Паншарпенинг при помощи 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

Мы видим, что картинка, производимая нашим алгоритмом выглядит как-будто она обработана фильтром, увеличивающим резкость. Кроме того, на ней легче различить оттенки, чем на той, которая является результатом работы OTB. Если копать глубже, то выяснится, что OTB сохранил диапазон значений исходного (очевидно, происходит нормализация данных на одном из этапов обработки), тогда как после использования нашего скрипта диапазон значений пикселей мультиспектрального снимка изменился кардинально. Это означает, что растр после паншарпенинга в OTB потенциально может быть использован для последующей радиометрической калибровки, а после паншарпенинга нашим скриптом - нет. Однако указанные различия в результатах работы нашего скрипта и OTB следует считать не недостатками, а особенностями, которые следует учитывать в работе.

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

Ссылки по теме

Паншарпенинг в QGIS с использованием Orfeo ToolBox

Паншарпенинг в ERDAS

Обработка и интерпретация данных Landsat 8 (OLI) средствами GRASS GIS 7 [раздел Паншарпенинг в GRASS]

Геопроцессинг с SEXTANTE для QGIS