Использование консольного DTclassifier для классификации растровых данных и анализа изменений: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Нет описания правки
Нет описания правки
Строка 123: Строка 123:
</syntaxhighlight>
</syntaxhighlight>


Здесь следует пояснить, что при создании модели программа использует только те полигоны обучающей выборки которые пересекаются с указанными --input_rasters (если точнее, для всех остальных объектов, значения спектральных каналов будут равны нулю). Количество каналов (bands) по которым строится модель равно сумме всех каналов сцен перечисленных в --input_rasters. Значения пикселей пересекающие полигоны (также можно использовать линии или точки) обучающей выборки (тренингов), на основе которых и строится модель классификации, можно сохранить в виде точечного слоя (опция  --save_points). Также можно сохранить и саму модель (опция --save_model). Для каждой поддиректории (path/row) модель будет уникальной, поскольку строится только на тех полигонах (из слоя cut1.shp) которые пересекаются со снимком.  
Здесь следует пояснить, что при создании модели программа использует только те полигоны обучающей выборки которые пересекаются с указанными --input_rasters (если точнее, для всех остальных объектов, значения спектральных каналов будут равны нулю). Полигоны взяты из общего слоя, но для построения модели важны только те которые пересекаются с экстентом указанных в командной строке сцен. Количество каналов (bands) по которым строится модель равно сумме всех каналов сцен перечисленных в --input_rasters. Значения пикселей пересекающие полигоны (также можно использовать линии или точки) обучающей выборки (тренингов), на основе которых и строится модель классификации, можно сохранить в виде точечного слоя (опция  --save_points). Также можно сохранить и саму модель (опция --save_model). Для каждой поддиректории (path/row) модель будет уникальной, поскольку строится только на тех полигонах (из слоя cut1.shp) которые пересекаются со снимком.  
Для обработки нескольких path/row можно создать .bat файл, в котором последовательно записать командные строки для обработки всех директорий
Для обработки нескольких path/row можно создать .bat файл, в котором последовательно записать командные строки для обработки всех директорий
или создать скрипт (например на python):
или создать скрипт (например на python):

Версия от 19:00, 5 декабря 2016

Эта страница является черновиком статьи.


Описание и примеры использования консольной версии DTclassifier для классификации растровых данных.

DT classifier - это простой в использовании и эффективный плагин для классификации растровых изображений. Декстопная версия расширения входит в дистрибутив NextGIS QGIS. Расширение использует метод деревьев решений, реализованный на основе библиотеки OpenCV. Подробнее об установке и работе декстопной версии расширения написано здесь.

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


Тестовые данные

Загрузить архив c данными, использовавшимися при подготовке статьи (4.7Гб).

Работа с расширением

Синтаксис расширения включает несколько обязательных параметров и ряд дополнительных, расширяющих функциональность приложения.

DTconsol pic2a.jpg
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).

Создание обучающей выборки (тренингов) по снимкам Landsat до (декабрь 2015) и после (январь 2016) вырубки

Распакуем и сохраним на диске тестовые данные. В примере используется расположение 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:

Runscript.jpg

После выполнения операции результаты будут сохранены в директорию d:\DT\result (генерализованые растры имеют суффикс _smooth.tif). Также будут сохранены векторные слои (точки) треннингов (<path/row>train_points.shp) и модели (<path/row>model.yaml).


Результат классификации:

Результат анализа изменений по снимкам Landsat слева - общий результат, справа примеры снимков до (декабрь 2015) и после (январь-февраль 2016) нарушения, и результат классификации


Пример выявления выборочной рубки



Пример 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:

Результат классификации. Сверху-вниз примеры выявления выборочных рубок на снимках Landsat до и после нарушения. Дата снимка - в верхнем левом углу.

Цветом показаны результаты классификации на основе объединенной модели (красным) и модели для каждой сцены (синим - модель из предыдущего примера). Следует отметить, что результат объединенной модели несколько хуже (меньшая площадь рубок попала в изменения), что скорее всего связано с варьированием яркости между ненормализованными снимками.



Пример 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)
}

Результат автоматической обработки:

Результат классификации - маска облаков (синим), справа примеры снимка с результатом классификации (справа внизу) и без (справа вверху)

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

Данное расширение, благодаря простоте и относительно высокой производительности, позволит развивать свои проекты в области ДЗЗ не только профессиональным разработчикам, но и исследователям, только начинающим использовать языки программирования в своей работе.