Использование консольного DTclassifier для классификации растровых данных и анализа изменений: различия между версиями
Igg (обсуждение | вклад) Нет описания правки |
Нет описания правки |
||
(не показаны 92 промежуточные версии 5 участников) | |||
Строка 2: | Строка 2: | ||
{{Аннотация|Описание и примеры использования консольной версии DTclassifier для классификации растровых данных.}} | {{Аннотация|Описание и примеры использования консольной версии DTclassifier для классификации растровых данных.}} | ||
'''DTclassifier''' - это простой в использовании и эффективный плагин для классификации растровых изображений. | |||
Декстопная версия расширения входит в дистрибутив | Декстопная версия расширения входит в дистрибутив [http://nextgis.ru/nextgis-qgis/ NextGIS QGIS]. | ||
Расширение использует метод деревьев решений реализованный на основе библиотеки [ | Расширение использует метод деревьев решений, реализованный на основе библиотеки [https://github.com/opencv/opencv/wiki OpenCV]. Основная информация об установке и работе содержится в [http://gis-lab.info/qa/dtclassifier.html отдельной статье]. Если вы еще не знакомы с DTclassifier, рекомендуется начать его изучение с неё. | ||
Плагин имеет удобный графический интерфейс, но для некоторых задач, например, | Плагин имеет удобный графический интерфейс, но для некоторых задач, например, для встраивания алгоритмов классификации в процесс обработки изображений или для анализа большого количества данных, удобнее использовать консольную версию инструмента. Про это и пойдет речь в этой статье. | ||
для встраивания алгоритмов классификации в | |||
консольную версию | |||
==== Тестовые данные ==== | ==== Тестовые данные ==== | ||
[ | [https://drive.google.com/file/d/0B-Sx9ICoBa2QZzlCR2tOZEgtWkE/view?usp=sharing Загрузить] архив c данными, использовавшимися при подготовке статьи (4.7Гб). | ||
=== Работа с расширением === | === Работа с расширением === | ||
Строка 49: | Строка 17: | ||
расширяющих функциональность приложения. | расширяющих функциональность приложения. | ||
[[Файл:DTconsol_pic2a.jpg|700px|thumb|center]] | |||
<pre>classifier</pre> | <pre>classifier.bat</pre> | ||
- запуск расширения | - запуск расширения | ||
'''Основные опции''' | |||
<pre>--input_rasters</pre> | <pre>--input_rasters</pre> | ||
- пути к растровым данным ('''ВАЖНО'''! необходимо указывать полные пути к файлам) разделенные пробелом | - пути к растровым данным ('''ВАЖНО'''! необходимо указывать полные пути к файлам), разделенные пробелом | ||
<pre>--presence </pre> | <pre>--presence </pre> | ||
- полные пути к векторным слоям объектов содержащие признак который нужно выделить при анализе, разделенные пробелом. Например для природных экосистем: облака, водные объекты, рубки или пожары. | - полные пути к векторным слоям объектов, содержащие признак, который нужно выделить при анализе, разделенные пробелом. Например, для природных экосистем: облака, водные объекты, рубки или пожары. | ||
<pre>--absence</pre> | <pre>--absence</pre> | ||
- полные пути к векторным слоям фоновых объектов , т.е. объектам от которых нужно отделить объекты содержащие признак. | - полные пути к векторным слоям фоновых объектов, т.е. объектам, от которых нужно отделить объекты, содержащие признак. | ||
<pre>--classify</pre> | <pre>--classify</pre> | ||
путь для сохранения результата классификации | путь для сохранения результата классификации | ||
'''Дополнительные опции''' | |||
<pre>--decision_tree</pre> | |||
- использует дерево решений (по умолчанию используется случаный лес, т.е. random forest). ПРИМЕЧАНИЕ: Опыт показывает, что использование Random Forest обычно существенно улучшает результат классификации по сравнению с деревом решений | |||
<pre>--discrete_classes</pre> | |||
- Применяется при использовании одиночного дерева (--use_decision_tree). В этом случае, если флаг установлен, исходные данные будут трактоваться как набор дискретных величин (классификационная модель). Если флаг не установлен, то используется регрессионная модель. | |||
<pre>--generalize</pre> | |||
- размер окна для генерализации (kernel size). «Cглаживание» результатов классификации с настраиваемым размером окна. При включенном сглаживании будет создан не один, а два растра: классифицированный и сглаженный (имя содержит суффикс «_smooth»). | |||
<pre>-- | <pre>--save_train_layer </pre> | ||
- сохраняет | - сохраняет слой тренингов, полученный из обучающих данных (полигонов, линий и тд.) в виде точечного шейп-файла с атрибутивной таблицей, включающей значения растров и классы обучающей выборки (1- presence,0 - absence) | ||
<pre>--save_model</pre> | <pre>--save_model</pre> | ||
- | - сохраняет модель (дерево решений или модель random forest) в файл (.yaml) | ||
<pre>--use_model</pre> | <pre>--use_model</pre> | ||
- | - использует существующую модель. При этом параметры --presence и --absence игнорируются. | ||
('''ВАЖНО'''! количество и порядок каналов должен быть таким же как при создании модели, иначе результат будет некорректным) | ('''ВАЖНО'''! количество и порядок каналов должен быть таким же как при создании модели, иначе результат будет некорректным) | ||
<pre>-- | <pre>--use_train_layer</pre> | ||
- использует существующий точечный файл (полученный при помощи опции --save_train_layer) для создания модели. При этом игнорируются опции --presence и --absence. Также используется для создания модели без последующей классификации. В этом случае достаточно просто игнорировать опции --classify и --input_rasters | |||
< | Также приведены примеры синтаксиса командной строки: | ||
Классификация. | |||
<syntaxhighlight lang="python"> | |||
classifier.bat --input_rasters rast1 [rast2, ...] --presence vect1.shp [vect2, ...] --absence vect1.shp [vect2, ...] --classify result.tiff | |||
</syntaxhighlight> | |||
Классификация c с использованием существующей модели. | |||
<syntaxhighlight lang="python"> | |||
classifier.bat --input_rasters rast1 [rast2, ...] --use_model model.yaml --classify result.tiff | |||
</syntaxhighlight> | |||
< | Создание модели (без классификации). | ||
<syntaxhighlight lang="python"> | |||
classifier.bat --input_rasters rast1 [rast2, ...] --presence vect1.shp [vect2, ...] --absence vect1.shp [vect2, ...] --save_model model.yaml | |||
</syntaxhighlight> | |||
Создание модели на основе точечного слоя. | |||
<syntaxhighlight lang="python"> | |||
classifier.bat --use_train_layer train_points.shp --save_model model.yaml | |||
</syntaxhighlight> | |||
< | Создание точечного слоя тренингов. | ||
<syntaxhighlight lang="python"> | |||
classifier.bat --input_rasters rast1 [rast2, ...] --presence vect1.shp [vect2, ...] --absence vect1.shp [vect2, ...] --save_train_layer train_points.shp | |||
</syntaxhighlight> | |||
'''TODO''': Проверить! Если слой тренингов выходит за границы растра, то часть точек за границами создается с ненулевыми значениями (+ errors при выполнении), которые приходится удалять вручную. | |||
=== Примеры использования === | === Примеры использования === | ||
==== Классификация сцен Landsat в batch режиме на основе общего слоя тестовых объектов ==== | ==== Пример 1. Классификация сцен Landsat в batch режиме на основе общего слоя тестовых объектов ==== | ||
Для выявления изменения лесного покрова в результате рубок в зимний период используется анализ изменений 2-х и более сцен Landsat. Если территория исследования занимает значительную площадь | Для выявления изменения лесного покрова в результате рубок в зимний период используется анализ изменений 2-х и более сцен Landsat. Если территория исследования занимает значительную площадь | ||
и включает в себя несколько path/row спутника, то удобнее использовать консольную версию для проведения классификации в batch режиме. В данном примере, сцены Landsat за период 2015-2016 гг. | и включает в себя несколько path/row спутника, то удобнее использовать консольную версию для проведения классификации в batch режиме. В данном примере, сцены Landsat за период 2015-2016 гг. | ||
разбиты по директориям и классификация проводится последовательно для каждого path/row. | разбиты по директориям (data/<path/row>) и классификация проводится последовательно для каждого path/row. | ||
Обучающая выборка (тренинги) состоит из двух шейп-файлов (сut1.shp - рубки, произошедшие за период наблюдения и nochange.shp - не изменившиеся территории). Тренинги находятся в отдельной директории (trainings). | |||
[[Файл:Pic1.jpg|700px|thumb|center|Создание обучающей выборки (тренингов) по снимкам Landsat до (декабрь 2015) и после (январь 2016) вырубки]] | |||
Распакуем и сохраним на диске тестовые данные. В примере используется расположение d:\DT. | |||
Снимки лежат в d:\DT\data, тренинги в d:\DT\trainings. Командная строка для классификации отдельного path/row: | |||
<syntaxhighlight lang="python"> | |||
classifier.bat --presence d:\DT\trainings\cut1.shp --absence d:\DT\trainings\nochange.shp --input_rasters d:\DT\data\113027\113027_16055.tif d:\DT\data\113027\113027_15020.tif --save_points --save_model d:\DT\data\113027\model.aml --classify d:\DT\data\113027\113027_dtclass.tif --generalize 3 | |||
</syntaxhighlight> | |||
Здесь следует пояснить, что при создании модели программа использует только те полигоны обучающей выборки которые пересекаются с указанными --input_rasters (если точнее, для всех остальных объектов, значения спектральных каналов будут равны нулю). Полигоны взяты из общего слоя, но для построения модели важны только те которые пересекаются с экстентом указанных в командной строке сцен. Количество каналов (bands) по которым строится модель равно сумме всех каналов сцен перечисленных в --input_rasters. Значения пикселей пересекающие полигоны (также можно использовать линии или точки) обучающей выборки (тренингов), на основе которых и строится модель классификации, можно сохранить в виде точечного слоя (опция --save_points). Также можно сохранить и саму модель (опция --save_model). Для каждой поддиректории (path/row) модель будет уникальной, поскольку строится только на тех полигонах (из слоя cut1.shp) которые пересекаются со снимком. | |||
Для обработки нескольких path/row можно создать .bat файл, в котором последовательно записать командные строки для обработки всех директорий | |||
или создать скрипт (например на python): | |||
Пример [https://github.com/IgorGlushkov/DT_classifier_utilities/blob/master/DT_batch_classify.py скрипта] для обработки нескольких сцен | |||
<syntaxhighlight lang="python"> | |||
# -*- coding: UTF-8 -*- | |||
import os | |||
import glob | |||
# path to data | |||
datadir='d:\DT\data' | |||
# path to training | |||
traindir='d:\DT\\trainings' | |||
# presence | |||
presence=os.path.join(traindir,'cut1.shp') | |||
# absence | |||
absence=os.path.join(traindir,'nochange.shp') | |||
#result | |||
result='d:\DT\\result' | |||
for subdir, dirs, files in os.walk(datadir): | |||
scenes=glob.glob(os.path.join(subdir, '*.tif')) | |||
if scenes == []: | |||
continue | |||
else: | |||
scenes=' '.join(scenes) | |||
output=os.path.join(result, str(os.path.basename(subdir))+'out.tif') | |||
model=os.path.join(result, str(os.path.basename(subdir))+'model.yaml') | |||
train_points=os.path.join(result, str(os.path.basename(subdir))+'train_points.shp') | |||
# command DTclassifier | |||
command= 'classifier.bat --input_rasters %s --save_train_layer %s --save_model %s --presence %s --absence %s --classify %s --generalize 3' % (scenes,train_points,model,presence,absence,output) | |||
#run command | |||
os.system(command) | |||
}</syntaxhighlight> | |||
Запустим скрипт в OSGeo4W shell: | |||
[[Файл:Runscript.jpg|700px|thumb|center]] | |||
После выполнения операции результаты будут сохранены в директорию d:\DT\result (генерализованые растры имеют суффикс _smooth.tif). Также будут сохранены векторные слои (точки) треннингов (<path/row>train_points.shp) и модели (<path/row>model.yaml). | |||
Результат классификации: | |||
[[Файл:DTconsol_pic4.jpg|700px|thumb|center|Результат анализа изменений по снимкам Landsat слева - общий результат, справа примеры снимков до (декабрь 2015) и после (январь-февраль 2016) нарушения, и результат классификации]] | |||
[[Файл:Pic2res.jpg|700px|thumb|center|Пример выявления выборочной рубки]] | |||
==== Пример 2. Создание объединенной модели классификации ==== | |||
ВАЖНО: Этот пример носит скорее демонстрационный характер, поскольку для создания реальной модели входящие данные должны быть [http://gis-lab.info/qa/regress-r.html нормализованы]. | |||
Предыдущий пример фактически повторяет функционал декстопной версии, а в данном примере используются две новые функции DTClassifier ----use_train_layer и --use_model, которые позволяют собирать тренировочные данные с отдельно выбранных сцен и на этой основе создавать модель, которую затем можно применить ко всему массиву данных. | |||
Структура входящих данных (число каналов или сцен) должна быть одинаковой для каждого блока данных, в отличие от предыдущего примера, где число снимков / каналов в поддиректории может варьировать, поскольку для каждого path/row создается отдельная модель для классификации. В нашем случае количество каналов для каждого path/row одинаковое (2 снимка по 7 каналов в каждой поддиректории), поэтому в данном примере мы используем сохраненные на предыдущем этапе точечные данные (--save_points) для создания общей для всех сцен модели и последовательно применим ее к каждой поддиректории с данными. | |||
Для этого создадим объединенный слой точек train_points.shp, просто объединив ..113027\out_train_points.shp и ..112027\out_train_points.shp в QGIS. Дополнительно можно убрать нулевые значения | |||
что выровняет баланс между числом тренировочных точек, содержащих признак, и точками с его отсутствием. | |||
[[Файл:DTconsol pic6.jpeg|700px|thumb|center|Точечный слой треннингов для созданий модели]] | |||
Для отдельной рубки слой выглядит так: | |||
[[Файл:Ex2 trainings1.jpg|700px|thumb|center|Точечный слой треннингов для отдельной рубки (пример)]] | |||
Сохраним объединенный слой точек в директорию d:\DT\model\points.shp. Затем создадим модель, используя командную строку: | |||
<syntaxhighlight lang="python"> | |||
classifier.bat --use_train_layer d:\DT\model\points.shp --save_model d:\DT\model\model_all.yaml | |||
</syntaxhighlight> | |||
И применим ее последовательно к данным, к примеру, используя скрипт (python): | |||
<syntaxhighlight lang="python"> | |||
# -*- coding: UTF-8 -*- | |||
import os | |||
import glob | |||
# path to data | |||
datadir='d:\DT\data' | |||
# path to model | |||
model='d:\DT\model\model_all.yaml' | |||
#result | |||
result='d:\DT\\result' | |||
for subdir, dirs, files in os.walk(datadir): | |||
scenes=glob.glob(os.path.join(subdir, '*.tif')) | |||
if scenes == []: | |||
continue | |||
else: | |||
scenes=' '.join(scenes) | |||
output=os.path.join(result, str(os.path.basename(subdir))+'out1.tif') | |||
#command DTclassifier | |||
command= 'classifier.bat --input_rasters %s --use_model %s --classify %s --generalize 3' % (scenes,model,output) | |||
#run command | |||
print command | |||
os.system(command) | |||
}</syntaxhighlight> | |||
После завершения, в директории d:\DT\result добавятся результаты классификации (<path/row>out1.tif и <path/row>out1_smooth.tif). | |||
Откроем результат в QGIS: | |||
[[Файл:Ex2 resmodel.jpg|700px|thumb|center|Результат классификации. Сверху-вниз примеры выявления выборочных рубок на снимках Landsat до и после нарушения. Дата снимка - в верхнем левом углу.]] | |||
Цветом показаны результаты классификации на основе объединенной модели (красным) и модели для каждой сцены (синим - модель из предыдущего примера). | |||
Следует отметить, что результат объединенной модели несколько хуже (меньшая площадь рубок попала в изменения), что скорее всего связано с варьированием яркости между ненормализованными снимками. | |||
==== Пример 3. Cоздание модели облачности на основе ''исходных'' данных Landsat ==== | |||
'''ВАЖНО''': Этот пример также носит скорее демонстрационный характер, поскольку для создания реальной модели каналы лучше сначала перевести в значения излучения на сенсоре ([https://grasswiki.osgeo.org/wiki/Atmospheric_correction toar reflectance]). | |||
Данный пример иллюстрирует возможность встраивания плагина в процесс обработки данных, например создание масок облачности (или любой другой классификации) для нескольких сцен, на основе общего слоя треннингов. | |||
В директории DT\cloud_data находится три "исходных" сцены Landsat (..tar.gz - архивы), а в DT\cloud_trainings обучающие слои облаков (cloud.shp) и необлачных участков (nocloud.shp). В данном примере слои получены случайной выборкой из классифицированного до этого (другими методами) растра, но также можно использовать любые другие способы создания тестовых объектов. | |||
Попробуем максимально автоматизировать процесс обработки и получения масок облачности для данных сцен. | |||
Для этого воспользуемся возможностями python, целиком скрипт доступен [https://github.com/IgorGlushkov/DT_classifier_utilities/blob/master/DT_cloud_model.py здесь]. | |||
Разобьем обработку на несколько этапов: | |||
'''Определение путей к данным и тренингам и импорт необходимых библиотек.''' | |||
<syntaxhighlight lang="python"> | |||
# -*- coding: UTF-8 -*- | |||
import os | |||
import glob | |||
import tarfile | |||
# path to data | |||
datadir='d:\DT\cloud_data' | |||
# path to training | |||
traindir='d:\DT\cloud_trainings' | |||
# presence | |||
presence=os.path.join(traindir,'cloud.shp') | |||
# absence | |||
absence=os.path.join(traindir,'nocloud.shp') | |||
#merge points | |||
merge_train_points='cloud_train_points.shp' | |||
#model | |||
model='d:\DT\model\model_cloud.yaml' | |||
#result | |||
result='d:\DT\\result' | |||
}</syntaxhighlight> | |||
'''Распаковка и удаление ненужных каналов.''' | |||
В данном примере я просто удалил панхроматические и BQA каналы, но конечно можно использовать и другие опции (и код), чтобы сохранить эти данные для других задач. | |||
<syntaxhighlight lang="python"> | |||
#unzip in separate folders | |||
scenes = glob.glob(os.path.join(datadir, '*gz')) | |||
for scene in scenes: | |||
try: | |||
a = tarfile.open(scene) | |||
scene_name = scene.split('.')[0] | |||
a.extractall(path=os.path.join(datadir,scene_name)) | |||
a.close() | |||
except IOError: | |||
False | |||
#remove BQA bands | |||
for subdir, dirs, files in os.walk(datadir): | |||
bands=glob.glob(os.path.join(subdir, '*BQA.TIF')) | |||
if bands == []: | |||
continue | |||
else: | |||
command= 'rm %s' % (bands[0]) | |||
print command | |||
os.system(command) | |||
#remove panchromatic bands | |||
for subdir, dirs, files in os.walk(datadir): | |||
bands=glob.glob(os.path.join(subdir, '*B8.TIF')) | |||
if bands == []: | |||
continue | |||
else: | |||
command= 'rm %s' % (bands[0]) | |||
print command | |||
os.system(command) | |||
}</syntaxhighlight> | |||
'''Создание общего точечного слоя тренингов.''' | |||
Здесь следует отметит, что создание обучающего слоя тренингов - ключевой момент, который определяет качество получаемой модели и классификации. Для данного примера важно, чтобы тренинги с облаками не располагались в области перекрытия снимков, поскольку в этом случае возникнут дублирующиеся точки с разными значениями класса объекта, но с одинаковыми значениями излучения в каналах снимка. | |||
Также следует учитывать, что DT classifier собирает данные со всех тренингов последовательно для каждого снимка, что приводит к появлению большого числа нулевых значений там, где тренинги выходят за границы снимка, поэтому при объединении задана функция "использовать только ненулевые значения каналов (в данном примере -where Band_1 > 0)". При использовании большего числа снимков необходимо будет учитывать и дату каждой отдельной сцены при создании тренингов, но в данном простом примере используются только неперекрывающиеся участки облачности. | |||
<syntaxhighlight lang="python"> | |||
#create train points from selected scenes | |||
for subdir, dirs, files in os.walk(datadir): | |||
bands=glob.glob(os.path.join(subdir, '*.TIF')) | |||
if bands == []: | |||
continue | |||
else: | |||
bands=' '.join(bands) | |||
train_points=os.path.join(traindir, str(os.path.basename(subdir))+'train_points.shp') | |||
command= 'classifier.bat --input_rasters %s --save_train_layer %s --presence %s --absence %s' % (bands,train_points,presence,absence) | |||
print command | |||
os.system(command) | |||
#merge trainings with null values removing using ogr2ogr | |||
commandlist=[] | |||
driver='ESRI Shapefile' | |||
trainings=glob.glob(os.path.join(traindir, '*train_points.shp')) | |||
for training in trainings: | |||
if commandlist == []: | |||
merge_train=os.path.join(traindir, merge_train_points) | |||
command='ogr2ogr -f \"%s\" %s %s -where \"Band_1 > 0\"' % (driver,merge_train,training) | |||
merge_name = merge_train_points.split('.')[0] | |||
command1='ogr2ogr -f \"%s\" -update -append %s %s -nln %s -where \"Band_1 > 0\"' % (driver,merge_train,training,merge_name) | |||
commandlist.append(command) | |||
commandlist.append(command1) | |||
else: | |||
merge_train=os.path.join(traindir, merge_train_points) | |||
#remove ext | |||
merge_name = merge_train_points.split('.')[0] | |||
command='ogr2ogr -f \"%s\" -update -append %s %s -nln %s -where \"Band_1 > 0\"' % (driver,merge_train,training,merge_name) | |||
commandlist.append(command) | |||
#run merge shp | |||
#f = open('%s\com.txt' % (traindir), 'w') | |||
for command in commandlist: | |||
#f.write(command+'\n') | |||
os.system(command) | |||
#f.close() | |||
}</syntaxhighlight> | |||
В результате - создан точечный слой d:\DT\cloud_trainings\cloud_train_points.shp. | |||
[[Файл:Ex3 trainigs fin.jpg|700px|thumb|center|Точечный слой треннингов на облачность]] | |||
'''Создание модели и классификация.''' | |||
Заключительный этап - создание модели из объединенного слоя тренингов и классификация. | |||
Результаты сохранены в директории ..DT\result в двух файлах для каждого снимка, <имя исходного архива>..out_cloud.tif и ..out_cloud_smooth.tif (генерализованный) | |||
<syntaxhighlight lang="python"> | |||
#create model | |||
if os.path.isfile(os.path.join(traindir,merge_train_points)): | |||
try: | |||
merge_train=os.path.join(traindir,merge_train_points) | |||
command='classifier.bat --use_train_layer %s --save_model %s' % (merge_train,model) | |||
os.system(command) | |||
except IOError: | |||
False | |||
#classify | |||
for subdir, dirs, files in os.walk(datadir): | |||
bands=glob.glob(os.path.join(subdir, '*.TIF')) | |||
if bands == []: | |||
continue | |||
else: | |||
bands=' '.join(bands) | |||
output=os.path.join(result, str(os.path.basename(subdir))+'out_cloud.tif') | |||
command= 'classifier.bat --input_rasters %s --use_model %s --classify %s --generalize 3' % (bands,model,output) | |||
print command | |||
os.system(command) | |||
}</syntaxhighlight> | |||
Результат автоматической обработки: | |||
[[Файл:Ex3 result fin.jpg|700px|thumb|center|Результат классификации - маска облаков (синим), справа примеры снимка с результатом классификации (справа внизу) и без (справа вверху)]] | |||
В качестве заключения можно сказать, что автоматизация задач обработки и классификации данных ДЗЗ развивается стремительными темпами, и хотя по-прежнему остается прерогативой крупных команд и организаций, | |||
методы полу-автоматической и автоматической обработки становятся доступны и для более широкого круга пользователей. | |||
Данное расширение, благодаря простоте и относительно высокой производительности, позволит развивать свои проекты в области ДЗЗ не только профессиональным разработчикам, но и исследователям, только начинающим использовать языки программирования в своей работе. |
Текущая версия от 19:59, 11 декабря 2016
Описание и примеры использования консольной версии DTclassifier для классификации растровых данных.
DTclassifier - это простой в использовании и эффективный плагин для классификации растровых изображений. Декстопная версия расширения входит в дистрибутив NextGIS QGIS. Расширение использует метод деревьев решений, реализованный на основе библиотеки OpenCV. Основная информация об установке и работе содержится в отдельной статье. Если вы еще не знакомы с DTclassifier, рекомендуется начать его изучение с неё.
Плагин имеет удобный графический интерфейс, но для некоторых задач, например, для встраивания алгоритмов классификации в процесс обработки изображений или для анализа большого количества данных, удобнее использовать консольную версию инструмента. Про это и пойдет речь в этой статье.
Тестовые данные
Загрузить архив c данными, использовавшимися при подготовке статьи (4.7Гб).
Работа с расширением
Синтаксис расширения включает несколько обязательных параметров и ряд дополнительных, расширяющих функциональность приложения.
classifier.bat
- запуск расширения
Основные опции
--input_rasters
- пути к растровым данным (ВАЖНО! необходимо указывать полные пути к файлам), разделенные пробелом
--presence
- полные пути к векторным слоям объектов, содержащие признак, который нужно выделить при анализе, разделенные пробелом. Например, для природных экосистем: облака, водные объекты, рубки или пожары.
--absence
- полные пути к векторным слоям фоновых объектов, т.е. объектам, от которых нужно отделить объекты, содержащие признак.
--classify
путь для сохранения результата классификации
Дополнительные опции
--decision_tree
- использует дерево решений (по умолчанию используется случаный лес, т.е. random forest). ПРИМЕЧАНИЕ: Опыт показывает, что использование Random Forest обычно существенно улучшает результат классификации по сравнению с деревом решений
--discrete_classes
- Применяется при использовании одиночного дерева (--use_decision_tree). В этом случае, если флаг установлен, исходные данные будут трактоваться как набор дискретных величин (классификационная модель). Если флаг не установлен, то используется регрессионная модель.
--generalize
- размер окна для генерализации (kernel size). «Cглаживание» результатов классификации с настраиваемым размером окна. При включенном сглаживании будет создан не один, а два растра: классифицированный и сглаженный (имя содержит суффикс «_smooth»).
--save_train_layer
- сохраняет слой тренингов, полученный из обучающих данных (полигонов, линий и тд.) в виде точечного шейп-файла с атрибутивной таблицей, включающей значения растров и классы обучающей выборки (1- presence,0 - absence)
--save_model
- сохраняет модель (дерево решений или модель random forest) в файл (.yaml)
--use_model
- использует существующую модель. При этом параметры --presence и --absence игнорируются. (ВАЖНО! количество и порядок каналов должен быть таким же как при создании модели, иначе результат будет некорректным)
--use_train_layer
- использует существующий точечный файл (полученный при помощи опции --save_train_layer) для создания модели. При этом игнорируются опции --presence и --absence. Также используется для создания модели без последующей классификации. В этом случае достаточно просто игнорировать опции --classify и --input_rasters
Также приведены примеры синтаксиса командной строки:
Классификация.
classifier.bat --input_rasters rast1 [rast2, ...] --presence vect1.shp [vect2, ...] --absence vect1.shp [vect2, ...] --classify result.tiff
Классификация c с использованием существующей модели.
classifier.bat --input_rasters rast1 [rast2, ...] --use_model model.yaml --classify result.tiff
Создание модели (без классификации).
classifier.bat --input_rasters rast1 [rast2, ...] --presence vect1.shp [vect2, ...] --absence vect1.shp [vect2, ...] --save_model model.yaml
Создание модели на основе точечного слоя.
classifier.bat --use_train_layer train_points.shp --save_model model.yaml
Создание точечного слоя тренингов.
classifier.bat --input_rasters rast1 [rast2, ...] --presence vect1.shp [vect2, ...] --absence vect1.shp [vect2, ...] --save_train_layer train_points.shp
TODO: Проверить! Если слой тренингов выходит за границы растра, то часть точек за границами создается с ненулевыми значениями (+ errors при выполнении), которые приходится удалять вручную.
Примеры использования
Пример 1. Классификация сцен Landsat в batch режиме на основе общего слоя тестовых объектов
Для выявления изменения лесного покрова в результате рубок в зимний период используется анализ изменений 2-х и более сцен Landsat. Если территория исследования занимает значительную площадь и включает в себя несколько path/row спутника, то удобнее использовать консольную версию для проведения классификации в batch режиме. В данном примере, сцены Landsat за период 2015-2016 гг. разбиты по директориям (data/<path/row>) и классификация проводится последовательно для каждого path/row. Обучающая выборка (тренинги) состоит из двух шейп-файлов (сut1.shp - рубки, произошедшие за период наблюдения и nochange.shp - не изменившиеся территории). Тренинги находятся в отдельной директории (trainings).
Распакуем и сохраним на диске тестовые данные. В примере используется расположение d:\DT. Снимки лежат в d:\DT\data, тренинги в d:\DT\trainings. Командная строка для классификации отдельного path/row:
classifier.bat --presence d:\DT\trainings\cut1.shp --absence d:\DT\trainings\nochange.shp --input_rasters d:\DT\data\113027\113027_16055.tif d:\DT\data\113027\113027_15020.tif --save_points --save_model d:\DT\data\113027\model.aml --classify d:\DT\data\113027\113027_dtclass.tif --generalize 3
Здесь следует пояснить, что при создании модели программа использует только те полигоны обучающей выборки которые пересекаются с указанными --input_rasters (если точнее, для всех остальных объектов, значения спектральных каналов будут равны нулю). Полигоны взяты из общего слоя, но для построения модели важны только те которые пересекаются с экстентом указанных в командной строке сцен. Количество каналов (bands) по которым строится модель равно сумме всех каналов сцен перечисленных в --input_rasters. Значения пикселей пересекающие полигоны (также можно использовать линии или точки) обучающей выборки (тренингов), на основе которых и строится модель классификации, можно сохранить в виде точечного слоя (опция --save_points). Также можно сохранить и саму модель (опция --save_model). Для каждой поддиректории (path/row) модель будет уникальной, поскольку строится только на тех полигонах (из слоя cut1.shp) которые пересекаются со снимком.
Для обработки нескольких path/row можно создать .bat файл, в котором последовательно записать командные строки для обработки всех директорий
или создать скрипт (например на python):
Пример скрипта для обработки нескольких сцен
# -*- coding: UTF-8 -*-
import os
import glob
# path to data
datadir='d:\DT\data'
# path to training
traindir='d:\DT\\trainings'
# presence
presence=os.path.join(traindir,'cut1.shp')
# absence
absence=os.path.join(traindir,'nochange.shp')
#result
result='d:\DT\\result'
for subdir, dirs, files in os.walk(datadir):
scenes=glob.glob(os.path.join(subdir, '*.tif'))
if scenes == []:
continue
else:
scenes=' '.join(scenes)
output=os.path.join(result, str(os.path.basename(subdir))+'out.tif')
model=os.path.join(result, str(os.path.basename(subdir))+'model.yaml')
train_points=os.path.join(result, str(os.path.basename(subdir))+'train_points.shp')
# command DTclassifier
command= 'classifier.bat --input_rasters %s --save_train_layer %s --save_model %s --presence %s --absence %s --classify %s --generalize 3' % (scenes,train_points,model,presence,absence,output)
#run command
os.system(command)
}
Запустим скрипт в OSGeo4W shell:
После выполнения операции результаты будут сохранены в директорию d:\DT\result (генерализованые растры имеют суффикс _smooth.tif). Также будут сохранены векторные слои (точки) треннингов (<path/row>train_points.shp) и модели (<path/row>model.yaml).
Результат классификации:
Пример 2. Создание объединенной модели классификации
ВАЖНО: Этот пример носит скорее демонстрационный характер, поскольку для создания реальной модели входящие данные должны быть нормализованы.
Предыдущий пример фактически повторяет функционал декстопной версии, а в данном примере используются две новые функции DTClassifier ----use_train_layer и --use_model, которые позволяют собирать тренировочные данные с отдельно выбранных сцен и на этой основе создавать модель, которую затем можно применить ко всему массиву данных.
Структура входящих данных (число каналов или сцен) должна быть одинаковой для каждого блока данных, в отличие от предыдущего примера, где число снимков / каналов в поддиректории может варьировать, поскольку для каждого path/row создается отдельная модель для классификации. В нашем случае количество каналов для каждого path/row одинаковое (2 снимка по 7 каналов в каждой поддиректории), поэтому в данном примере мы используем сохраненные на предыдущем этапе точечные данные (--save_points) для создания общей для всех сцен модели и последовательно применим ее к каждой поддиректории с данными.
Для этого создадим объединенный слой точек train_points.shp, просто объединив ..113027\out_train_points.shp и ..112027\out_train_points.shp в QGIS. Дополнительно можно убрать нулевые значения что выровняет баланс между числом тренировочных точек, содержащих признак, и точками с его отсутствием.
Для отдельной рубки слой выглядит так:
Сохраним объединенный слой точек в директорию d:\DT\model\points.shp. Затем создадим модель, используя командную строку:
classifier.bat --use_train_layer d:\DT\model\points.shp --save_model d:\DT\model\model_all.yaml
И применим ее последовательно к данным, к примеру, используя скрипт (python):
# -*- coding: UTF-8 -*-
import os
import glob
# path to data
datadir='d:\DT\data'
# path to model
model='d:\DT\model\model_all.yaml'
#result
result='d:\DT\\result'
for subdir, dirs, files in os.walk(datadir):
scenes=glob.glob(os.path.join(subdir, '*.tif'))
if scenes == []:
continue
else:
scenes=' '.join(scenes)
output=os.path.join(result, str(os.path.basename(subdir))+'out1.tif')
#command DTclassifier
command= 'classifier.bat --input_rasters %s --use_model %s --classify %s --generalize 3' % (scenes,model,output)
#run command
print command
os.system(command)
}
После завершения, в директории d:\DT\result добавятся результаты классификации (<path/row>out1.tif и <path/row>out1_smooth.tif).
Откроем результат в QGIS:
Цветом показаны результаты классификации на основе объединенной модели (красным) и модели для каждой сцены (синим - модель из предыдущего примера). Следует отметить, что результат объединенной модели несколько хуже (меньшая площадь рубок попала в изменения), что скорее всего связано с варьированием яркости между ненормализованными снимками.
Пример 3. Cоздание модели облачности на основе исходных данных Landsat
ВАЖНО: Этот пример также носит скорее демонстрационный характер, поскольку для создания реальной модели каналы лучше сначала перевести в значения излучения на сенсоре (toar reflectance).
Данный пример иллюстрирует возможность встраивания плагина в процесс обработки данных, например создание масок облачности (или любой другой классификации) для нескольких сцен, на основе общего слоя треннингов. В директории DT\cloud_data находится три "исходных" сцены Landsat (..tar.gz - архивы), а в DT\cloud_trainings обучающие слои облаков (cloud.shp) и необлачных участков (nocloud.shp). В данном примере слои получены случайной выборкой из классифицированного до этого (другими методами) растра, но также можно использовать любые другие способы создания тестовых объектов. Попробуем максимально автоматизировать процесс обработки и получения масок облачности для данных сцен.
Для этого воспользуемся возможностями python, целиком скрипт доступен здесь.
Разобьем обработку на несколько этапов:
Определение путей к данным и тренингам и импорт необходимых библиотек.
# -*- coding: UTF-8 -*-
import os
import glob
import tarfile
# path to data
datadir='d:\DT\cloud_data'
# path to training
traindir='d:\DT\cloud_trainings'
# presence
presence=os.path.join(traindir,'cloud.shp')
# absence
absence=os.path.join(traindir,'nocloud.shp')
#merge points
merge_train_points='cloud_train_points.shp'
#model
model='d:\DT\model\model_cloud.yaml'
#result
result='d:\DT\\result'
}
Распаковка и удаление ненужных каналов.
В данном примере я просто удалил панхроматические и BQA каналы, но конечно можно использовать и другие опции (и код), чтобы сохранить эти данные для других задач.
#unzip in separate folders
scenes = glob.glob(os.path.join(datadir, '*gz'))
for scene in scenes:
try:
a = tarfile.open(scene)
scene_name = scene.split('.')[0]
a.extractall(path=os.path.join(datadir,scene_name))
a.close()
except IOError:
False
#remove BQA bands
for subdir, dirs, files in os.walk(datadir):
bands=glob.glob(os.path.join(subdir, '*BQA.TIF'))
if bands == []:
continue
else:
command= 'rm %s' % (bands[0])
print command
os.system(command)
#remove panchromatic bands
for subdir, dirs, files in os.walk(datadir):
bands=glob.glob(os.path.join(subdir, '*B8.TIF'))
if bands == []:
continue
else:
command= 'rm %s' % (bands[0])
print command
os.system(command)
}
Создание общего точечного слоя тренингов.
Здесь следует отметит, что создание обучающего слоя тренингов - ключевой момент, который определяет качество получаемой модели и классификации. Для данного примера важно, чтобы тренинги с облаками не располагались в области перекрытия снимков, поскольку в этом случае возникнут дублирующиеся точки с разными значениями класса объекта, но с одинаковыми значениями излучения в каналах снимка.
Также следует учитывать, что DT classifier собирает данные со всех тренингов последовательно для каждого снимка, что приводит к появлению большого числа нулевых значений там, где тренинги выходят за границы снимка, поэтому при объединении задана функция "использовать только ненулевые значения каналов (в данном примере -where Band_1 > 0)". При использовании большего числа снимков необходимо будет учитывать и дату каждой отдельной сцены при создании тренингов, но в данном простом примере используются только неперекрывающиеся участки облачности.
#create train points from selected scenes
for subdir, dirs, files in os.walk(datadir):
bands=glob.glob(os.path.join(subdir, '*.TIF'))
if bands == []:
continue
else:
bands=' '.join(bands)
train_points=os.path.join(traindir, str(os.path.basename(subdir))+'train_points.shp')
command= 'classifier.bat --input_rasters %s --save_train_layer %s --presence %s --absence %s' % (bands,train_points,presence,absence)
print command
os.system(command)
#merge trainings with null values removing using ogr2ogr
commandlist=[]
driver='ESRI Shapefile'
trainings=glob.glob(os.path.join(traindir, '*train_points.shp'))
for training in trainings:
if commandlist == []:
merge_train=os.path.join(traindir, merge_train_points)
command='ogr2ogr -f \"%s\" %s %s -where \"Band_1 > 0\"' % (driver,merge_train,training)
merge_name = merge_train_points.split('.')[0]
command1='ogr2ogr -f \"%s\" -update -append %s %s -nln %s -where \"Band_1 > 0\"' % (driver,merge_train,training,merge_name)
commandlist.append(command)
commandlist.append(command1)
else:
merge_train=os.path.join(traindir, merge_train_points)
#remove ext
merge_name = merge_train_points.split('.')[0]
command='ogr2ogr -f \"%s\" -update -append %s %s -nln %s -where \"Band_1 > 0\"' % (driver,merge_train,training,merge_name)
commandlist.append(command)
#run merge shp
#f = open('%s\com.txt' % (traindir), 'w')
for command in commandlist:
#f.write(command+'\n')
os.system(command)
#f.close()
}
В результате - создан точечный слой d:\DT\cloud_trainings\cloud_train_points.shp.
Создание модели и классификация.
Заключительный этап - создание модели из объединенного слоя тренингов и классификация. Результаты сохранены в директории ..DT\result в двух файлах для каждого снимка, <имя исходного архива>..out_cloud.tif и ..out_cloud_smooth.tif (генерализованный)
#create model
if os.path.isfile(os.path.join(traindir,merge_train_points)):
try:
merge_train=os.path.join(traindir,merge_train_points)
command='classifier.bat --use_train_layer %s --save_model %s' % (merge_train,model)
os.system(command)
except IOError:
False
#classify
for subdir, dirs, files in os.walk(datadir):
bands=glob.glob(os.path.join(subdir, '*.TIF'))
if bands == []:
continue
else:
bands=' '.join(bands)
output=os.path.join(result, str(os.path.basename(subdir))+'out_cloud.tif')
command= 'classifier.bat --input_rasters %s --use_model %s --classify %s --generalize 3' % (bands,model,output)
print command
os.system(command)
}
Результат автоматической обработки:
В качестве заключения можно сказать, что автоматизация задач обработки и классификации данных ДЗЗ развивается стремительными темпами, и хотя по-прежнему остается прерогативой крупных команд и организаций, методы полу-автоматической и автоматической обработки становятся доступны и для более широкого круга пользователей.
Данное расширение, благодаря простоте и относительно высокой производительности, позволит развивать свои проекты в области ДЗЗ не только профессиональным разработчикам, но и исследователям, только начинающим использовать языки программирования в своей работе.