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

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

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