Паншарпенинг в QGIS с использованием Orfeo ToolBox
В данной статье рассказывается о том как произвести паншарпенинг в QGIS посредством инструментария Orfeo ToolBox.
Введение
При работе с данными дистанционного зондирования Земли зачастую приходится решать задачу улучшения снимков с низким разрешением за счёт снимков с высоким. В QGIS нет встроенного инструмента паншарпенинга, но есть возможность подключения инструментов Orfeo ToolBox (OTB), где они есть.
Немного теории
На удивление сложно найти определение для слова "паншарпенинг". Попробуем сформулировать его своими словами, взяв за основу описание из книги Шовергердт Р.А. Дистанционное зондирование. Модели и методы обработки изображений - М.: Техносфера, 2010. (с. 416):
Паншарпенинг (от анг. Panchromatic sharpening) - это процесс объединения изображений в пространственной области основная задача которого заключается в передаче высокочастотного содержания изображения с высоким разрешением (обычно панхроматического) изображению с низким разрешением (обычно мультиспектральному).
Одним из наиболее простых алгоритмов паншарпенинга является модуляция высоких частот. В этом алгоритме для получения улучшенного изображения в канале PXS панхроматический снимок PAN попиксельно умножается на изображение с низким разрешением XS, после чего полученный результат нормируется на низкочастотную компоненту панхроматического снимка PANsmooth (панхроматический снимок, обработанный сглаживающим фильтром с окном, соответствующим размеру пикселя изображения с низким разрешением):
где i и j - индексы пикселей.
Паншарпенинг в OTB и QGIS
С программной точки зрения процесс паншарпенинга делится на две стадии: 1) подготовительную, во время которой разрешение и экстент мультиспектрального растра приводится в соответствие с экстентом и разрешением панхроматического растра (это необходимо для проведения операций растровой алгебры, вовлекающих оба растра); 2) непосредственно паншарпенинг.
Использование утилиты OTB
OTB по своей сути является набором консольных утилит и приложений. За паншарпенинг отвечает приложение otbcli_BundleToPerfectSensor. Оно совмещает в себе обе стадии, описанные выше, и запускается следующей консольной командой:
otbcli_BundleToPerfectSensor -inp pan_image -inxs xs_image -out output_image
где pan_image, xs_image и output_image соответсвенно - пути к панхроматическому, мультиспектральному и результирующему растрам. О дополнительных опциях читайте в документации.
Запуск паншарпенинга из QGIS
В QGIS инструменты OTB подключаются через модуль Processing. После установки OTB на компьютер зайдите в QGIS в настройки модуля Processing (Processing -> Options -> Providers) и активируйте интеграцию с OTB:
Теперь в Processing Toolbox у вас появится набор инструментов Orfeo Toolbox. К сожалению, в инструментарий входят только консольные утилиты, но не приложения (то есть модуль otbcli_BundleToPerfectSensor не доступен из QGIS). Поэтому не спешите сразу запускать утилиты паншарпенинга (Orfeo Toolbox -> Geomentry -> Pansharpening), так как здесь требуется, чтобы все входные растры имели идентичный экстент и разрешение, то есть выполняется только вторая стадия. Чтобы выполнить первую - необходимо воспользоваться модулем Superimpose sensor (Orfeo Toolbox -> Geomentry -> Superimpose sensor).
Модуль Superimpose sensor имеет следующие параметры:
- Reference input - базовое изображение, к параметрам которого будет приведено обрабатываемое изображение (здесь выбираем панхроматическое).
- The image to reproject - изображение, чьи параметры будут приведены к параметрам базового изображения (здесь выбираем мультиспектральное).
- Default elevation - средняя высота изображения над эллипсоидом.
- Spacing of the deformation field - интервал расстояния поля деформации (coarser deformation field).
- Interpolation - тип интерполяции.
- Available RAM - ограничение использования оперативной памяти.
- Output image - результирующий растр.
После трансформации мультиспектрального растра можно переходить непосредственно к паншарпенингу. Модуль Pansharpening(rsc) имеет следующие параметры:
- Input PAN image - панхроматическое изображение высокого разрешения.
- Input XS image - мультиспектральное изображение, приведённое к параметрам панхроматического.
- Algorithm - алгоритм для паншарпенинга (не изменяется, для использования других алгоритмов (LMVM и Bayesian) используйте соответсвующие модули Pansharpening(bayes) и Pansharpening(lvmv))
- Available RAM - ограничение использования оперативной памяти.
- Output image - результирующий растр.
Давайте посмотрим на результат паншарпенинга при помощи указанных модулей.
Результат выглядит довольно неплохо, и что немаловажно, был получен очень бытро - OTB использует многопоточную обработку и задействует все ядра вашего процессора.
Создание модели
Согласитесь, гораздо удобнее было бы обойтись без промежуточных операций и растров. Этого можно добиться путём создания соответствующей модели для Processing Toolbox. В качестве входных параметров модели следует использовать панхроматический и мультиспектральный растры. Последний сначала подготавливается при помощи Superimpose sensor, а затем обрабатывается при помощи Pansharpening(rsc).
Интерфейс модели будет выглядеть следующим образом:
- High resolution panchrome raster - панхроматическое изображение высокого разрешения.
- Low resolution raster - мультиспектральное изображение низкого разрешения.
- pansharpened OTB - результирующий растр.
Код модели
Ниже вы найдёте код, описывающий модель, приведённую выше. Вы можете сохранить его в виде текстового файла с расширением .model (например, Pansharp_OTB.model) и сохранить в папку ".../.qgis2/processing/models/". После этого она появится у вас в списке моделей Processing toolbox (необходим рестарт QGIS).
{
"values": {
"inputs": {
"RASTERLAYER_HIGHRESOLUTIONPANCHROMATICRASTER": {
"values": {
"pos": {
"values": {
"y": 55.0,
"x": 121.0
},
"class": "point"
},
"param": {
"values": {
"isAdvanced": false,
"name": "RASTERLAYER_HIGHRESOLUTIONPANCHROMATICRASTER",
"showSublayersDialog": true,
"value": null,
"exported": null,
"hidden": false,
"optional": false,
"description": "High_resolution_panchromatic_raster"
},
"class": "processing.core.parameters.ParameterRaster"
}
},
"class": "processing.modeler.ModelerAlgorithm.ModelerParameter"
},
"RASTERLAYER_LOWRESOLUTIONMULTISPECTRALRASTER": {
"values": {
"pos": {
"values": {
"y": 55.0,
"x": 342.0
},
"class": "point"
},
"param": {
"values": {
"isAdvanced": false,
"name": "RASTERLAYER_LOWRESOLUTIONMULTISPECTRALRASTER",
"showSublayersDialog": true,
"value": null,
"exported": null,
"hidden": false,
"optional": false,
"description": "Low_resolution_multispectral_raster"
},
"class": "processing.core.parameters.ParameterRaster"
}
},
"class": "processing.modeler.ModelerAlgorithm.ModelerParameter"
}
},
"helpContent": {},
"group": "Raster processing",
"name": "OTB Pan-sharpening",
"algs": {
"OTBPANSHARPENINGRCS_1": {
"values": {
"name": "OTBPANSHARPENINGRCS_1",
"paramsFolded": true,
"outputs": {
"-out": {
"values": {
"description": "Panshrpened_OTB",
"pos": {
"values": {
"y": 311.0,
"x": 556.0
},
"class": "point"
}
},
"class": "processing.modeler.ModelerAlgorithm.ModelerOutput"
}
},
"pos": {
"values": {
"y": 266.0,
"x": 356.0
},
"class": "point"
},
"outputsFolded": true,
"dependencies": [],
"params": {
"-inp": {
"values": {
"name": "RASTERLAYER_HIGHRESOLUTIONPANCHROMATICRASTER"
},
"class": "processing.modeler.ModelerAlgorithm.ValueFromInput"
},
"-inxs": {
"values": {
"alg": "OTBSUPERIMPOSESENSOR_1",
"output": "-out"
},
"class": "processing.modeler.ModelerAlgorithm.ValueFromOutput"
},
"-ram": 2048.0,
"-method": 0
},
"active": true,
"consoleName": "otb:pansharpeningrcs",
"description": "Pansharpening (rcs)"
},
"class": "processing.modeler.ModelerAlgorithm.Algorithm"
},
"OTBSUPERIMPOSESENSOR_1": {
"values": {
"name": "OTBSUPERIMPOSESENSOR_1",
"paramsFolded": true,
"outputs": {},
"pos": {
"values": {
"y": 182.0,
"x": 133.0
},
"class": "point"
},
"outputsFolded": true,
"dependencies": [],
"params": {
"-interpolator": 0,
"-inm": {
"values": {
"name": "RASTERLAYER_LOWRESOLUTIONMULTISPECTRALRASTER"
},
"class": "processing.modeler.ModelerAlgorithm.ValueFromInput"
},
"-ram": 2048.0,
"-interpolator.bco.radius": 2.0,
"-lms": 4.0,
"-inr": {
"values": {
"name": "RASTERLAYER_HIGHRESOLUTIONPANCHROMATICRASTER"
},
"class": "processing.modeler.ModelerAlgorithm.ValueFromInput"
},
"-elev.default": 0.0
},
"active": true,
"consoleName": "otb:superimposesensor",
"description": "Superimpose sensor"
},
"class": "processing.modeler.ModelerAlgorithm.Algorithm"
}
}
},
"class": "processing.modeler.ModelerAlgorithm.ModelerAlgorithm"
}
Ссылки по теме
Обработка и интерпретация данных Landsat 8 (OLI) средствами GRASS GIS 7