{{Статья|Черновик}}
{{Аннотация|В данной статье рассказывается о том как произвести паншарпенинг в программной среде R.

[via [http://ssrebelious.blogspot.ru/2015/02/pan-sharpening-using-r.html Misanthrope's Thoughts]]}}

==Введение==
Одной из наиболее рапространённых задач при работе со спутниковыми снимками является паншарпенинг. С этой задачёй успешно справляются и проприетарные и открытые программы для работы с ДДЗЗ. Но в некоторых случаях вам может захотеться иметь больший контроль над проводиммыми операциями или, быть может использовать алгоритм, ещё не реализованный ни в одном ПО. В таком случае на помощь придёт один из языков программирования, например, [https://ru.wikipedia.org/wiki/R_%28%D1%8F%D0%B7%D1%8B%D0%BA_%D0%BF%D1%80%D0%BE%D0%B3%D1%80%D0%B0%D0%BC%D0%BC%D0%B8%D1%80%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D1%8F%29 R].

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

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

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

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

[[Файл:Formula OTB pansharpening.png|center]]

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

==Реализация==
Для работы с растрами, имеющими географическую привязку в ''R'' есть замечательный пакет [http://cran.r-project.org/web/packages/raster/index.html ''raster'']. На основании функций, доступных в этом пакете мы и разработаем наш скрипт. В использовании он будет чрезвычайно прост:
<syntaxhighlight lang="rsplus">
library(raster)

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

Начнём создавать необходимые функции. Первая функция - самая простая. Она всего-лишь реализовывает формулу, представленную выше.
<syntaxhighlight lang="rsplus">
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]
}
</syntaxhighlight>

Вторая функция будет вычислять низкочастотную составляющую растра с высоким разрешением.
<syntaxhighlight lang="rsplus">
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'

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

Теперь создадим функцию, объединяющую две предыдущие и выдающую растр с улучшенным разрешением.
<syntaxhighlight lang="rsplus">
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) # итоговый растр
}
</syntaxhighlight>

Ну, и наконец, последняя функция, носящая чисто утилитарный характер (позволяет использовать минимум параметров при сохранении растров) - обёртка к функции сохранения растра из пакета 'raster'. Идея состоит в том, что мы почти всегда будем хотеть сохранить результат в формате GeoTIFF, и при этом иметь целочисленный тип данных. Дело в том, что в большинстве случаев значения пикселей мультиспектрального растра - целые числа, а после паншарпенинга они становятся числами с плавающей запятой. Десятичные значения пикселя не несут никакой полезной информации, а занимают лишние сотни мегабайт места на диске. Использование целых чисел в качестве типа данных растра экономит более 0,5 Гб дискового пространства при паншарпенинге сцен WorldView-2.
<syntaxhighlight lang="rsplus">
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)
}
</syntaxhighlight>

==Результаты и обсуждение==
[[Файл:Initial panchrome.png|thumb|center|x700px|Пахроматический растр с пространственным разрешением 0,5 м]]
[[Файл:Initial multi.png|thumb|center|x700px|Мультиспектральный растр с пространственным разрешением 2,0 м]]
[[Файл:R pansharp result.png|thumb|center|x700px|Мультиспектральный растр после паншарпенинга в R]]
Результат можно сравнить с результатом паншарпенинга в Orfeo TollBox:
[[Файл:OTB pansharp result.png|thumb|center|x700px|Мультиспектральный растр после паншарпенинга при помощи OTB]]
К сожалению, в raster не реализована многопоточность, поэтому производительность операций по обработке больших растров оставляет желать лучшего.

==Скрипт для Processing Tollbox в QGIS==
==Заключение==

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

2015-04-18T21:46:22Z

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

Файл:R pansharp result.png

2015-04-18T21:45:06Z

2015-04-04T12:51:25Z

SSSRebelious: растр после паншарпенинга в OTB

растр после паншарпенинга в OTB

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

2015-04-04T12:50:14Z

SSSRebelious: /* Паншарпенинг в OTB и QGIS */

Файл:Initial multi.png

2015-04-04T12:27:04Z

SSSRebelious: изначальный панхроматический растр

изначальный панхроматический растр