Просмотр кода страницы Распознавание сгоревших территорий с помощью деревьев решений и OpenCV

{{Статья|Опубликована|burnedarea-opencv}}
{{Аннотация|Использование OpenCV для работы со снимками MODIS и их классификации.}}

Лето 2010 года наглядно показало, что для нас с вами очень актуальна проблема природных пожаров. Для устранения их последствий необходимо знать, на каких именно территориях они произошли, знать их точное пространственное расположение. Один из перспективных подходов к решению этой задачи — классификация территорий с помощью [http://www.basegroup.ru/library/analysis/tree/description/ деревьев решений] по [http://gis-lab.info/qa/modislandprod.html данным MODIS – Land]. Для эффективной обработки растровых данных и их классификации хорошо подходит библиотека компьютерного зрения [http://opencv.willowgarage.com/ OpenCV] (Open Source Computer Vision). В ней эффективно реализована необходимая нам функциональность, библиотека бесплатная, свободная (лицензия BSD), кросс-платформенная и является фактически стандартом в компьютерном зрении — количество её скачиваний превысило 2 миллиона.

[[Файл:burnedarea-opencv-01.png|center]]

В этой статье мы по шагам рассмотрим как можно использовать библиотеку OpenCV для распознавания сгоревших территорий. Для примера мы будем использовать данные MODIS на начало и конец лета 2010 года представленными двумя продуктами:

* [https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/surface_reflectance/8_day_l3_global_500m/mod09a1 MOD09A1] — 8-дневные композиты, содержащие 7 каналов отражающей способности земной поверхности (Surface Reflectance Bands 1–7);
* [https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/burned_area/monthly_l3_global_500m/mcd45a1 MCD45A1] — информация о сгоревшей территории с указанием даты пожаров для каждого пикселя местности.

Оба этих продукта имеют разрешение 500 метров на один пиксель местности. В качестве тренировочной выборки мы возьмём данные по республике Марий-Эл, а в качестве тестовой — данные по Нижегородской области.

=== Импорт данных MODIS ===

Прежде всего нужно перевести данные MODIS из формата HDF в GeoTIFF — про это подробно написано в статье «[http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]». После этого мы сможем загрузить данные в OpenCV как обычное изображение в формате GeoTIFF. Поскольку изображение — это набор яркостей пикселей, то в OpenCV он будет хранится как матрица. Для этого используется класс Mat — главный класс для работы с матрицами и изображениями в OpenCV:

<syntaxhighlight lang="cpp">
Mat modisBand = imread("sur_refl_b01.tif", CV_LOAD_IMAGE_UNCHANGED);
//читаем изображение из файла "sur_refl_b01.tif" и записываем его в матрицу modisBand
</syntaxhighlight>

Подгрузив аналогично 7 каналов MODIS на начало лета и 7 каналов на конец лета, мы можем собрать все данные в одну 14-канальную матрицу, то есть прямоугольную матрицу элементами которой являются вектора из 14 элементов:

<syntaxhighlight lang="cpp">
vector<Mat> allBands;
allBands.push_back(modisBand);
... //подгружаем по очереди все каналы и добавляем их в allBands
Mat data;
merge(allBands, data); //объединяем все каналы в одну матрицу data
</syntaxhighlight>

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

<syntaxhighlight lang="cpp">
Mat burnDate = imread(burnDateFilename, CV_LOAD_IMAGE_UNCHANGED);
//элемент матрицы burnDate содержат информацию о дате пожара, произошедшего в соответствующем пикселе местности
Mat unburnedArea = (burnDate == 0);
//элемент матрицы unburnedArea будет равен 255, если пиксель не сгорел (дата пожара равна нулю), или 0 в противном случае.
uint16_t minBurnDay = 153; //2 июня
uint16_t maxBurnDay = 248; //5 сентября   
Mat burnedArea = (burnDate >= minBurnDay) & (burnDate <= maxBurnDay);
//элемент матрицы burnedArea будет равен 255, если в нём был пожар в указанном диапазоне дат, или 0 в противном случае
</syntaxhighlight>

=== Обучение деревьев решений ===

В качестве классификатора можно использовать простой и популярный подход — дерево решений. Но более точные результаты можно попробовать получить, если использовать [http://ru.wikipedia.org/wiki/Random_forest ансамбль деревьев решений (Random Forest)], который обобщает народную мудрость «одна голова хорошо, а две лучше». Сейчас данные загружены в виде 14-канальной матрицы. Для тренировки деревьев решений нам нужна 1-канальная матрица, в которой каждая строка — это один элемент тренировочной выборки, а в столбцах хранятся числовые значения его признаков (соответствующих каналов). То есть нам нужно преобразовать исходную матрицу в 1-канальную матрицу из 14 столбцов, в которой строка соответствует данным MODIS для одного пикселя местности:

<syntaxhighlight lang="cpp">
Mat trainData = data.reshape(1, data.total());
//сделать из data 1-канальную матрицу с числом строк, равным числу пикселей в data
Mat responses = burnedArea.reshape(1, burnedArea.total());
//аналогично поступаем с метками классов
</syntaxhighlight>

Теперь можем запустить тренировку ансамбля деревьев решений:

<syntaxhighlight lang="cpp">
CvRTrees classifier;
classifier.train(trainData, CV_ROW_SAMPLE, responses);
//СV_ROW_SAMPLE указывает, что элементы выборки для классификации хранятся по строкам
</syntaxhighlight>

Совершенно аналогично мы можем натренировать обычное дерево решений для сравнения:

<syntaxhighlight lang="cpp">
CvDTree decisionTree;
decisionTree.train(trainData, CV_ROW_SAMPLE, responses);
</syntaxhighlight>

=== Распознавание ===

Следующий шаг — классификация каждого пикселя местности тестовых данных с помощью натренированного ансамбля деревьев решений:

<syntaxhighlight lang="cpp">
int rowsCount = testData.rows; 
//количество строк в матрице testData, равное количеству классифицируемых пикселей местности
Mat prediction = Mat(rowsCount, 1, CV_32SC1); 
//создаём матрицу размера rowsCount x 1 и типом элемента int
for(int i=0; i<rowsCount; i++)
{
  Mat sample = testData.row(i); 
  //извлекаем i-ую строчку из матрицы testData - значения 14 каналов MODIS для соответствующего пикселя местности
  prediction.at<int>(i, 0) = classifier.predict(sample); 
  //записываем результат классификации в элемент (i, 0) матрицы prediction
}
</syntaxhighlight>

=== Пост-обработка ===

Деревья решений классифицируют каждый пиксель независимо друг от друга. Это приводит к тому, что хотя все соседи пикселя классифицированы как сгоревшие, он сам может может быть классифицирован как несгоревший — для реального распространения пожара это маловероятно. Для учёта результатов классификации соседних пикселей можно использовать операции [http://ru.wikipedia.org/wiki/%D0%9C%D0%B0%D1%82%D0%B5%D0%BC%D0%B0%D1%82%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B0%D1%8F_%D0%BC%D0%BE%D1%80%D1%84%D0%BE%D0%BB%D0%BE%D0%B3%D0%B8%D1%8F математической морфологии]:

<syntaxhighlight lang="cpp">
dilate(predictedBurnedArea, predictedBurnedArea, Mat(), Point(-1, -1), 5); //применяем наращивание 5 раз. Используется структурирующий элемент по умолчанию - квадратный элемент 3x3 erode(predictedBurnedArea, predictedBurnedArea, Mat(), Point(-1, -1), 7); //применяем эрозию 7 раз dilate(predictedBurnedArea, predictedBurnedArea, Mat(), Point(-1, -1), 3); //применяем наращивание 3 раза
</syntaxhighlight>

Результаты обработки (белым цветом показаны сгоревшие территории):

{|
|[[Файл:burnedarea-opencv-02.png|center]]
|[[Файл:burnedarea-opencv-03.png|center]]
|-
| до морфологии
| после морфологии
|}

=== Валидация ===

Для оценки качества классификации воспользуемся стандартным подходом — подсчитаем [http://gis-lab.info/qa/error-matrix.html матрицу ошибок (confusion matrix)]. Элемент <tt>(i, j)</tt> этой матрицы равен количеству объектов класса <tt>i</tt>, классифицированных как класс <tt>j</tt>. Матрицу ошибок можно посчитать, имея матрицы с реальными значениями классов (<tt>groundTruth</tt>) и с предсказанными значениями (<tt>prediction</tt>):

<syntaxhighlight lang="cpp">
int classesCount = 2; //количество распознаваемых классов
Mat confusionMatrix(classesCount, classesCount, CV_32SC1); 
//создаём матрицу размера classesCount x classesCount типа int
for(int i=0; i<classesCount; i++)
{  
  for(int j=0; j<classesCount; j++)
  {
    confusionMatrix.at<int>(i, j) = countNonZero((groundTruth == i) & (prediction == j)); 
    //подсчитываем количество элементов класса i, классифицированных как j    
  }
}
</syntaxhighlight>

=== Визуализация ===

Каналам RGB соответствуют каналы 1, 3 и 4 спутника MODIS. Для создания снимка нужно объединить эти каналы в одно изображение. Для повышения качества полученного изображения с ним нужно сделать несложные преобразования, описанные в статье «[http://gis-lab.info/docs/modis_true_color.pdf Creating Reprojected True Color MODIS Images: A Tutorial]». После этого отрисуем контуры сгоревших территорий поверх полученного изображения:

<syntaxhighlight lang="cpp">
vector<vector<Point>> contours; 
//массив для найденных контуров 
findContours(burnedArea, contours, CV_RETR_LIST, CV_CHAIN_APPROX_NONE); 
//найти контура по маске сгоревших территорий
drawContours(trueColorImage, contours, -1, Scalar(0, 0, 255), 2); 
//нарисовать все контуры красной линией толщины 2
namedWindow("burned area contours", CV_WINDOW_NORMAL); 
//создать окно c именем "burned area contours"
imshow("burned area contours", trueColorImage); 
//отобразить изображение в созданном окне
waitKey();
//отображать изображение, пока пользователь не нажмёт "any key"
</syntaxhighlight>

=== Экспорт результатов ===

Для экспорта результатов в GeoTIFF сначала сохраним полученное изображение в формате TIFF:

<syntaxhighlight lang="cpp">
imwrite("burnedArea.tif", trueColorImage);
</syntaxhighlight>

Теперь сохранённое изображение нужно связать с географическими координатами. Для этого есть несколько способов, в том числе есть и хорошие решения — кросс-платформенные, бесплатные и с открытым исходным кодом. Например, можно использовать библиотеку GDAL — про это подробно написано в статье «[http://gis-lab.info/qa/georef-gdal.html Использование GDAL для привязки растровых материалов]».

=== Результаты экспериментов ===

Для повышения точности классификации помимо описанных шагов лучше добавить ещё несколько дополнительных фильтров (например, фильтрация данных по качеству и тренировка классификатора только на надёжных данных). Полную версию алгоритма с данными и исходным кодом можно скачать [http://gis-lab.info/other/burnedarea-opencv.zip здесь]. Теперь мы можем визуально оценить результаты работы алгоритма:

[[Файл:burnedarea-opencv-04.png|center]]

Количественную оценку дают построенные матрицы ошибок:

Матрица ошибок для одного обычного дерева решений

{|class="wikitable"
|rowspan="2" colspan="2" |
|colspan="2" | Предсказанный класс территории
|-
|Несгоревшая
|Сгоревшая
|-
|rowspan="2" | Реальный класс территории
|Несгоревшая 	
|109423
|17387
|-
|Сгоревшая
|60
|3372
|}

Матрица ошибок для ансамбля деревьев решений

{|class="wikitable"
|rowspan="2" colspan="2" |
|colspan="2" | Предсказанный класс территории
|-
|Несгоревшая
|Сгоревшая
|-
|rowspan="2" | Реальный класс территории
|Несгоревшая 	
|120367
|6443
|-
|Сгоревшая
|172
|3260
|}

Рассчитаем ошибки [http://gis-lab.info/qa/error-matrix.html омиссии и комиссии]:

{|class="wikitable"
|rowspan="2" colspan="2" |
|colspan="2" | Классификатор
|-
|Дерево решений
|Ансамбль деревьев решений
|-
|rowspan="2" | Ошибка
|омиссии 	
|0.02
|0.05
|-
|комиссии
|0.84
|0.66
|}

Используя ансамбль деревьев решений, мы нашли 95% случившихся пожаров, но при этом 66% территорий, классифицированных как сгоревшие, относятся к несгоревшим по данным MCD45, которые в свою очередь отличаются большой ошибкой омиссии («[http://gis-lab.info/qa/mcd45.html Данные о сгоревших площадях MCD45: описание и получение]»).

Не менее важной характеристикой алгоритма является время его работы. На Intel Core i7 960 @ 3.20GHz с использованием только одного потока время тренировки ансамбля деревьев решений составило 0.85 секунды (54000 пикселей), время обработки тестовой территории — 0.13 секунды (130 000 пикселей).

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

=== Ссылки по теме ===

* [http://opencv.willowgarage.com/ Библиотека компьютерного зрения OpenCV]
* [http://gis-lab.info/qa/mcd45.html Данные о сгоревших площадях MCD45: описание и получение]
* [http://gis-lab.info/qa/firms.html Мониторинг пожаров на природных территориях с помощью сервиса FIRMS]