Обсудить в форуме Комментариев 7Редактировать в вики
Использование OpenCV для работы со снимками MODIS и их классификации.
Лето 2010 года наглядно показало, что для нас с вами очень актуальна проблема природных пожаров. Для устранения их последствий необходимо знать, на каких именно территориях они произошли, знать их точное пространственное расположение. Один из перспективных подходов к решению этой задачи — классификация территорий с помощью деревьев решений по данным MODIS – Land. Для эффективной обработки растровых данных и их классификации хорошо подходит библиотека компьютерного зрения OpenCV (Open Source Computer Vision). В ней эффективно реализована необходимая нам функциональность, библиотека бесплатная, свободная (лицензия BSD), кросс-платформенная и является фактически стандартом в компьютерном зрении — количество её скачиваний превысило 2 миллиона.
В этой статье мы по шагам рассмотрим как можно использовать библиотеку OpenCV для распознавания сгоревших территорий. Для примера мы будем использовать данные MODIS на начало и конец лета 2010 года представленными двумя продуктами:
Оба этих продукта имеют разрешение 500 метров на один пиксель местности. В качестве тренировочной выборки мы возьмём данные по республике Марий-Эл, а в качестве тестовой — данные по Нижегородской области.
Содержание |
Прежде всего нужно перевести данные MODIS из формата HDF в GeoTIFF — про это подробно написано в статье «Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT». После этого мы сможем загрузить данные в OpenCV как обычное изображение в формате GeoTIFF. Поскольку изображение — это набор яркостей пикселей, то в OpenCV он будет хранится как матрица. Для этого используется класс Mat — главный класс для работы с матрицами и изображениями в OpenCV:
Mat modisBand = imread("sur_refl_b01.tif", CV_LOAD_IMAGE_UNCHANGED); //читаем изображение из файла "sur_refl_b01.tif" и записываем его в матрицу modisBand
Подгрузив аналогично 7 каналов MODIS на начало лета и 7 каналов на конец лета, мы можем собрать все данные в одну 14-канальную матрицу, то есть прямоугольную матрицу элементами которой являются вектора из 14 элементов:
vector<Mat> allBands; allBands.push_back(modisBand); ... //подгружаем по очереди все каналы и добавляем их в allBands Mat data; merge(allBands, data); //объединяем все каналы в одну матрицу data
Для обучения деревьев решений и валидации полученных результатов нам необходима информация о том, какие же территории действительно сгорели за это лето. С помощью OpenCV мы легко сможем получить эту информацию из продукта MCD45A1, используя логические операции между матрицами и скалярами, которые выполняются попиксельно:
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 в противном случае
В качестве классификатора можно использовать простой и популярный подход — дерево решений. Но более точные результаты можно попробовать получить, если использовать ансамбль деревьев решений (Random Forest), который обобщает народную мудрость «одна голова хорошо, а две лучше». Сейчас данные загружены в виде 14-канальной матрицы. Для тренировки деревьев решений нам нужна 1-канальная матрица, в которой каждая строка — это один элемент тренировочной выборки, а в столбцах хранятся числовые значения его признаков (соответствующих каналов). То есть нам нужно преобразовать исходную матрицу в 1-канальную матрицу из 14 столбцов, в которой строка соответствует данным MODIS для одного пикселя местности:
Mat trainData = data.reshape(1, data.total()); //сделать из data 1-канальную матрицу с числом строк, равным числу пикселей в data Mat responses = burnedArea.reshape(1, burnedArea.total()); //аналогично поступаем с метками классов
Теперь можем запустить тренировку ансамбля деревьев решений:
CvRTrees classifier; classifier.train(trainData, CV_ROW_SAMPLE, responses); //СV_ROW_SAMPLE указывает, что элементы выборки для классификации хранятся по строкам
Совершенно аналогично мы можем натренировать обычное дерево решений для сравнения:
CvDTree decisionTree; decisionTree.train(trainData, CV_ROW_SAMPLE, responses);
Следующий шаг — классификация каждого пикселя местности тестовых данных с помощью натренированного ансамбля деревьев решений:
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 }
Деревья решений классифицируют каждый пиксель независимо друг от друга. Это приводит к тому, что хотя все соседи пикселя классифицированы как сгоревшие, он сам может может быть классифицирован как несгоревший — для реального распространения пожара это маловероятно. Для учёта результатов классификации соседних пикселей можно использовать операции математической морфологии:
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 раза
Результаты обработки (белым цветом показаны сгоревшие территории):
до морфологии | после морфологии |
Для оценки качества классификации воспользуемся стандартным подходом — подсчитаем матрицу ошибок (confusion matrix). Элемент (i, j) этой матрицы равен количеству объектов класса i, классифицированных как класс j. Матрицу ошибок можно посчитать, имея матрицы с реальными значениями классов (groundTruth) и с предсказанными значениями (prediction):
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 } }
Каналам RGB соответствуют каналы 1, 3 и 4 спутника MODIS. Для создания снимка нужно объединить эти каналы в одно изображение. Для повышения качества полученного изображения с ним нужно сделать несложные преобразования, описанные в статье «Creating Reprojected True Color MODIS Images: A Tutorial». После этого отрисуем контуры сгоревших территорий поверх полученного изображения:
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"
Для экспорта результатов в GeoTIFF сначала сохраним полученное изображение в формате TIFF:
imwrite("burnedArea.tif", trueColorImage);
Теперь сохранённое изображение нужно связать с географическими координатами. Для этого есть несколько способов, в том числе есть и хорошие решения — кросс-платформенные, бесплатные и с открытым исходным кодом. Например, можно использовать библиотеку GDAL — про это подробно написано в статье «Использование GDAL для привязки растровых материалов».
Для повышения точности классификации помимо описанных шагов лучше добавить ещё несколько дополнительных фильтров (например, фильтрация данных по качеству и тренировка классификатора только на надёжных данных). Полную версию алгоритма с данными и исходным кодом можно скачать здесь. Теперь мы можем визуально оценить результаты работы алгоритма:
Количественную оценку дают построенные матрицы ошибок:
Матрица ошибок для одного обычного дерева решений
Предсказанный класс территории | |||
Несгоревшая | Сгоревшая | ||
Реальный класс территории | Несгоревшая | 109423 | 17387 |
Сгоревшая | 60 | 3372 |
Матрица ошибок для ансамбля деревьев решений
Предсказанный класс территории | |||
Несгоревшая | Сгоревшая | ||
Реальный класс территории | Несгоревшая | 120367 | 6443 |
Сгоревшая | 172 | 3260 |
Рассчитаем ошибки омиссии и комиссии:
Классификатор | |||
Дерево решений | Ансамбль деревьев решений | ||
Ошибка | омиссии | 0.02 | 0.05 |
комиссии | 0.84 | 0.66 |
Используя ансамбль деревьев решений, мы нашли 95% случившихся пожаров, но при этом 66% территорий, классифицированных как сгоревшие, относятся к несгоревшим по данным MCD45, которые в свою очередь отличаются большой ошибкой омиссии («Данные о сгоревших площадях MCD45: описание и получение»).
Не менее важной характеристикой алгоритма является время его работы. На Intel Core i7 960 @ 3.20GHz с использованием только одного потока время тренировки ансамбля деревьев решений составило 0.85 секунды (54000 пикселей), время обработки тестовой территории — 0.13 секунды (130 000 пикселей).
Таким образом, используя стандартные средства OpenCV, мы смогли легко получить простой и быстрый алгоритм для распознавания сгоревших территорий.
Обсудить в форуме Комментариев 7Редактировать в вики
Последнее обновление: 2014-05-14 22:48
Дата создания: 31.03.2011
Автор(ы): Илья Лысенков (инженер по программному обеспечению, Itseez)
© GIS-Lab и авторы, 2002-2021. При использовании материалов сайта, ссылка на GIS-Lab и авторов обязательна. Содержание материалов - ответственность авторов. (подробнее).