Эта статья является продолжением цикла статей посвященных развитию стартапа "Arrow". Ребята из моей команды тоже не отстают и те, кого больше интересует бизнес-сторона вопроса можете почитать "Старт проекта и гибкость как залог успеха: путь команды ARROW", а те кто больше по разработке оставайтесь здесь. Первая техническая статья на примере нашего проекта освещает проблему построения Веб-ГИС приложения и сервисов. Логика приложения, в том числе, будет реализовывать инструменты получения спутниковых снимков и их обработки по специальным алгоритмам, что позволит решать задачи мониторинга территорий - дешифрировать, классифицировать и кластеризировать объекты и явления местности на основе технологий машинного обучения.
В начале небольшой "ликбез", буквально на три строчки и максимально доступно. Дело в том, ребятушки, что наши глазки видят весьма интересно, на самом деле они чувствуют внешние раздражители - электромагнитные волны. Наш мозг интерпретирует эти раздражения в зрительные образы, которым обучается с детства.
Одним из главных источников этих волн являются звезды и соответственно наше Солнце. Свет идущий от него - это целый спектр волн разной длины. Падая на м��стные предметы большая его часть поглощается, но некоторые волны отражаются и попадают к нам в глаза или фотографу на матрицу фотоаппарата, или на сенсор космического спутника. От травы отражается волна длиной, соответствующей зеленой части спектра, от помидора - красной и т.д. Но сами волны бесцветны, все что в них полезное, так это их длина и, зависящая от этой длины, спектральная яркость. Цвет объектам придает наш, раздражающийся этой волной через глаза, мозг, а в съемочной аппаратуре - специальный светофильтр и фотоматрица. Фильтр пропускает волны с длиной соответствующей красной части спектра R(red) или зеленой G(green) или синей B(blue). Волны воздействуют на фотодиоды матрицы. Это воздействие, пройдя еще некоторую обработку, преобразуется в цифровые значения спектральных яркостей.
Таким образом, получаются трехканальные (R, G, B) растровые наборы данных, где каждая ячейка матрицы содержит значения спектральных яркостей - это цветосинтезированное изображение в видимом диапазоне длин волн.
Весь фокус в том, что мы с вами, ввиду ограниченности возможностей наших природных сенсоров, не в состоянии чувствовать большую часть электромагнитного спектра, теряя много информации об окружающей нас действительности.
Если бы мы смогли зафиксировать отраженное излучение в других частях электромагнитного спектра, то узнали бы о территории гораздо больше. И мы можем, точнее это могут сенсоры космических аппаратов, имеющие в своем распоряжении не три, а большее число каналов, например у спутника Lansat8 их одинадцать.
Каналы космического спутника Lansat8:
Канал 1 (Coastal/Aerosol) – голубой и фиолетовый. Сильно рассеивается пылью и частицами атмосферы, используется для обнаружения дыма и тумана. Обеспечивает оценку «здоровья» биосферы океана по степени яркости.
Каналы 2, 3, 4 (Blue, Green, Red) – видимые цвета спектра (это то что видят наши глазки).
Канал 5 (Near Infrared, NIR) – невидимый ближний ИК канал. Обеспечивает оценку состояния растительности суши по степени яркости, что обусловлено сильным отражением здоровой растительности.
Каналы 6, 7 (Shortwave Infrared, SWIR-1/2) – инфракрасные каналы. Излучение сильно поглощается водой. Данные используются для разл��чения видов растительности и почвы, облаков, снега и льда.
Канал 8 (Panchromatic) обеспечивает максимальное пространственное разрешение 15 м/пиксел.
Канал 9 (Cirrus) обеспечивает различимость перистых и кучевых облаков.
Каналы 10, 11 (Thermal Infrared, TIR-1/2) – каналы дальнего ИК (теплового) диапазона, дающие информацию о температуре поверхности Земли.
Комбинируя эти каналы определенным образом человек может видеть, анализировать и классифицировать то, что до этого момента было скрыто от его глаз.
В зависимости от типа исследуемых по снимкам объектов используются различные комбинации спектральных каналов, обеспечивающие максимальную различимость объектов по отношению к фону. В ряде случаев комбинации каналов применяются последовательно с целью поэтапного отделения объектов и их свойств. Стандартная комбинация в «естественных цветах», рассмотренная ранее, подразумевает размещение каналов R, G, B в соответствующих позициях (слотах) многоканального растра, рассматриваемого на экране компьютера (визуальное дешифрирование), либо обрабатываемого машинными методами классификации (автоматизированное дешифрирование). При визуальном дешифрирова��ии на три слота растра могут помещаться любые спектральные каналы (особенно информативны ближний и средний ИК диапазоны), либо что-то еще – например, каналы разновременных снимков, результат расчета вегетационного индекса и т.п. При автоматизированном дешифрировании число каналов в многоканальном растре может быть больше трех.
Для исследований объектов на поверхности Земли наиболее востребованы каналы с номерами 1-7 с разрешением 30 м/пикс, общее число комбинаций которых составляет 343 (все мы их конечно же рассматривать не будем). Основные комбинации и их свойства представлены ниже.
На данном этапе работа реализована в Jupyter notebook, для работы с растровыми наборами используются три библиотеки Python: GDAL, numpy и matplotlib. Сами снимки (Landsat8) были получены с этого ресурса. У снимков перед анализом было повышено пространственное разрешение с помощью паншарпенинга до 15 м/пикс.
#Зададим рабочуют дирректорию
workdir = "./Lansat-8_2021_05_16"
# Импорт библиотек
from osgeo import gdal
from osgeo import ogr
# Для работы с матрицами используем библиотеку numpy
import numpy
# Для визуализации растра используем библиотеку matplotlib
import matplotlib.pyplot as pltПервым делом создадим композит в естественных цветах (тот самый про который я говорил в самом начале статьи). В нашем распоряжении три одноканальных растра (band_2, band_3, band_4) соответствующих отраженному сигналу в синем, зеленом и красном видимых диапазонах длин волн соответственно.
# Исходные растры
b2_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B2.TIF")
b3_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B3.TIF")
b4_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B4.TIF")
# Читаем значения ячеек в канале 2 как массив
b_data = b2_ds.GetRasterBand(1).ReadAsArray()
# Посмотрим на матрицу. Срез c 500 по 1500 строку, с 500 по 505 столбец
b_data[3500:4500, 1500:2505]На выходе мы увидим ту самую матрицу значений спектральных яркостей о которой говорилось ранее. Исходные растры 16-битные, это видно по значениям матрицы, максимальное из которых может быть 65536 с этим будут связаны некоторые проблемы о решении которых чуть ниже.
array([[7144, 7116, 7098, ..., 7297, 7334, 7321],
[7119, 7143, 7126, ..., 7278, 7309, 7317],
[7124, 7219, 7172, ..., 7247, 7275, 7288],
...,
[7430, 7390, 7374, ..., 6988, 7026, 7025],
[7427, 7402, 7362, ..., 7007, 7023, 7038],
[7399, 7396, 7372, ..., 7091, 7057, 7058]], dtype=uint16)Следующим шагом нужно создать пустой геопривязанный набор растровых данных в формате GeoTIFF, задать ему количество каналов которое хотим поместить в него, систему координат и проекцию в которых он будет существовать и его привязку.
# Создаем новый датасет
# 1. Создаем драйвер
driver = gdal.GetDriverByName('GTiff')
# 2. Создаем набор данных Composit_rgb.TIF, копируя основные метаданные из любого одиночного канала
new_dataset = driver.Create(workdir + r".\Lansat-8_2021_05_16\Composit_rgb.TIF", b2_ds.RasterXSize, b2_ds.RasterYSize, 3, b2_ds.GetRasterBand(1).DataType)
# 3. Также берем из исходника информацию о СК и привязке
new_dataset.SetProjection(b2_ds.GetProjection())
new_dataset.SetGeoTransform(b2_ds.GetGeoTransform())
# ИЛИ Установим свою трансформацию
#new_dataset.SetGeoTransform((452845.0, 60.0, 0.0, -1999245.0, 0.0, -60.0))
# 4. Также берем из исходников массивы данных
new_dataset.GetRasterBand(1).WriteArray(b2_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(2).WriteArray(b3_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(3).WriteArray(b4_ds.GetRasterBand(1).ReadAsArray())
# 5. удаляя созданную переменную мы заставляем GDAL записать датасет в файл Composit_rgb.TIF на диск
# Чистим оперативную память
del new_dataset
# ПРОСМОТР ЧЕРЕЗ ArcGis или QGISПосмотрим на то что получилось прямо в Jupyter (хотя можно для этого использовать интерфейс любой ГИС, здесь чуть ниже я часто буду представлять результат полученный с помощью ArcMap 10.8). Для этого нужно три канала растра объединить в один массив с помощью numpy и, так как растр 16-битный, а библиотека matplotlib работает с 8-битными, нормализовать полученный массив.
# Нормализация
def normalize(input_band):
min_value, max_value = input_band.min()*1.0, input_band.max()*1.0
return ((input_band*1.0 - min_value*1.0)/(max_value*1.0 - min_value*1.0))
# вызываем функцию для всех трех каналов
comp = gdal.Open(workdir + r".\Lansat-8_2021_05_16\Composit_rgb.TIF")
comp_red_band_array = comp.GetRasterBand(1).ReadAsArray()
comp_green_band_array = comp.GetRasterBand(2).ReadAsArray()
comp_blue_band_array = comp.GetRasterBand(3).ReadAsArray()
comp_red_band_array_norm = normalize(comp_red_band_array)
comp_green_band_array_norm = normalize(comp_green_band_array)
comp_blue_band_array_norm = normalize(comp_blue_band_array)
comp_rgb_norm = numpy.dstack([comp_red_band_array_norm, comp_green_band_array_norm, comp_blue_band_array_norm])
plt.rcParams['figure.figsize'] = (8, 8) # размер экстента
plt.close()
plt.imshow(comp_rgb_norm)
А так этот снимок можно оформить в ГИС, ну да ладно. Теперь же переместимся в мир невиданного и попробуем получить композиты на основе комбинаций других каналов, а заодно и интерпретировать результат. Интересно? Поехали.
Синтез каналов 5-4-3 (NIR, Red, Green)
# Исходные растры
nir_band_array = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B5.TIF")
red_band_arra = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B4.TIF")
green_band_array = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B3.TIF")
# Создаем новый датасет
# 1. Создаем драйвер
driver = gdal.GetDriverByName('GTiff')
# 2. Создаем набор данных, копируя основные метаданные из любого одиночного канала
new_dataset = driver.Create(workdir + r".\Lansat-8_2021_05_16\Composit_5_4_3.tif", nir_band_array.RasterXSize, nir_band_array.RasterYSize, 3, nir_band_array.GetRasterBand(1).DataType)
# 3. Также берем из исходника информацию о СК и привязке
new_dataset.SetProjection(nir_band_array.GetProjection())
new_dataset.SetGeoTransform(nir_band_array.GetGeoTransform())
# 4. Также берем из исходников массивы данных
new_dataset.GetRasterBand(1).WriteArray(nir_band_array.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(2).WriteArray(red_band_arra.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(3).WriteArray(green_band_array.GetRasterBand(1).ReadAsArray())
# 5. удаляя созданную переменную мы заставляем GDAL записать датасет в файл Composit_rgb.TIF на диск
# Чистим оперативную память
del new_dataset
# ПРОСМОТР ЧЕРЕЗ ArcGis или QGIS
Видели когда-нибудь такое буйство красок? :) Возрите же ! Ха-ха-ха !!!
Применяется для изучения состояния растительного покрова, мониторинга дренажа и почвенной мозаики, изучения сельскохозяйственных культур. Насыщенные оттенки красного являются индикаторами здоровой и (или) широколиственной растительности, в то время как более светлые оттенки характеризуют травянистую растительность или редколесья/кустарники. Растительность в этой комбинации имеет оттенки красного, городская застройка – зелено-голубые цвета, почва – от темно до светло коричневого или серого, лед, снег и облака – белые или светло голубые. Хвойные леса по сравнению с лиственными имеют более темно-красную или даже коричневую окраску. ТБО, как вариант могут также светиться голубым.
Это только начало, двигаемся дальше.
Синтез каналов 7-5-3 (SWIR, NIR, Green)
# Исходные растры
swir_band_array = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B7.TIF")
nir_band_arra = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B5.TIF")
green_band_array = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B3.TIF")
# Создаем новый датасет
# 1. Создаем драйвер
driver = gdal.GetDriverByName('GTiff')
# 2. Создаем набор данных, копируя основные метаданные из любого одиночного канала
new_dataset = driver.Create(workdir + r".\Lansat-8_2021_05_16\Composit_7_5_3.tif", swir_band_array.RasterXSize, swir_band_array.RasterYSize, 3, swir_band_array.GetRasterBand(1).DataType)
# 3. Также берем из исходника информацию о СК и привязке
new_dataset.SetProjection(swir_band_array.GetProjection())
new_dataset.SetGeoTransform(swir_band_array.GetGeoTransform())
# 4. Также берем из исходников массивы данных
new_dataset.GetRasterBand(1).WriteArray(swir_band_array.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(2).WriteArray(nir_band_arra.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(3).WriteArray(green_band_array.GetRasterBand(1).ReadAsArray())
# 5. удаляя созданную переменную мы заставляем GDAL записать датасет в файл Composit_rgb.TIF на диск
# Чистим оперативную память
del new_dataset
# ПРОСМОТР ЧЕРЕЗ ArcGis или QGIS
Ну как вам? Данная комбинация весьма полезна при анализе пустынь, может быть использована для изучения сельскохозяйственных земель и водно-болотных угодий. Пройденные пожарами территории выглядят ярко красными. Городская застройка отображается в оттенках розово-фиолетового. Здоровая растительность выглядит ярко зеленой, травянистые сообщества – зелеными, ярко розовые участки детектируют открытую почву, коричневые и оранжевые тона характерны для разреженной растительности. Сухостойная растительность выглядит оранжевой, вода – голубой, фиолетовой. Песок, открытая почва и минералы могут быть представлены большим числом цветов и оттенков.
Ну я думаю на счет кода вы поняли он однотипен меняем только каналы участвующие в синтезе, поэтому далее еще несколько примеров просто с интерпретацией результата.
Синтез каналов 5-6-2 (NIR, SWIR, Blue)
Добавление SWIR канала обеспечивает различимость возраста растительности. Здоровая растительность отображается в оттенках красного, коричневого, оранжевого и зеленого. Почвы могут выглядеть зелеными или коричневыми, урбанизированные территории – белесыми, серыми и зелено-голубыми, ярко голубой цвет может детектировать недавно вырубленные территории, а красноватые – восстановление растительности или разреженную растительность. Чистая, глубокая вода будет выглядеть темно синей (почти черной), для мелководья или высокого содержания взвесей в цвете преобладают более светлые синие оттенки.
Синтез каналов 7-6-4 (SWIR-2, SWIR-1, Red)
Применяется для мониторинга пожаров, так как тепловые аномалии выглядят красноватыми или желтыми. Также хорошо выделяются затопленные территории. Они имеют темно-синий и почти черный цвет. Поглощение излучения в среднем ИК диапазоне водой позволяет четко выделять береговую линию и водные объекты на снимке. Растительность отображается в оттенках темно и светло зеленого, урбанизированные территории выглядят белыми, зелено-голубыми и малиновыми, почвы, песок и минералы могут иметь много разных цветов.
Синтез каналов 6-5-4 (SWIR, NIR, Red)
Удобен для изучения растительного покрова и широко используется для анализа состояния лесных сообществ. Комбинация дает очень много информации и цветовых контрастов. Здоровая растительность выглядит ярко зеленой, а почвы – розовато-лиловыми. В отличие от синтеза 7-4-2, включающего канал SWIR2 и позволяющего изучать геологические процессы, эта комбинация дает возможность лучше различать и анализировать сельскохозяйственные угодья.
Ну и так далее, можно играть и выискивать различные результаты до бесконечности.
А что если нам нужно не просто классифицировать территорию, а еще и производить некие расчеты с применением полученного результата, т.е. работать с количественными данными? Тут нам на помощь прийдет "Алгебра растров". Например, нужно рассчитать нормализованный индекс растительности, давайте попробуем.
Нормализованный разностный индекс растительности (NDVI)
Это стандартизированный индекс, показывающий наличие и состояние растительности (относительную биомассу). Оценку необходимо производить в вегетационный период. Для расчета применяется следующая формула.
где NIR - значения пикселов из ближнего и Red - значения пикселов из красного канала.
Создадим семиканальный композит, почему, да потому что NDWI это не единственный спектральный индекс, который можно рассчитать, их много и в общей сложности их можно классифицировать на несколько следующих категорий.
Индексы растительности и почвы;
Индексы воды;
Геологические индексы;
Ландшафтные индексы.
# Создадим 7-канальный композит
# Исходные растры
b1_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B1.TIF")
b2_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B2.TIF")
b3_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B3.TIF")
b4_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B4.TIF")
b5_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B5.TIF")
b6_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B6.TIF")
b7_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\LC08_L1TP_188022_20211101_20211109_01_T1_B7.TIF")
# Создаем новый датасет
# 1. Создаем драйвер
driver = gdal.GetDriverByName('GTiff')
# 2. Создаем набор данных Composit_rgb.TIF, копируя основные метаданные из любого одиночного канала
new_dataset = driver.Create(workdir + r".\Lansat-8_2021_05_16\Composit_bands_1_7.TIF", b1_ds.RasterXSize, b1_ds.RasterYSize, 7, b1_ds.GetRasterBand(1).DataType)
# 3. Также берем из исходника информацию о СК и привязке
new_dataset.SetProjection(b1_ds.GetProjection())
new_dataset.SetGeoTransform(b1_ds.GetGeoTransform())
# 4. Также берем из исходников массивы данных
new_dataset.GetRasterBand(1).WriteArray(b1_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(2).WriteArray(b2_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(3).WriteArray(b3_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(4).WriteArray(b4_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(5).WriteArray(b5_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(6).WriteArray(b6_ds.GetRasterBand(1).ReadAsArray())
new_dataset.GetRasterBand(7).WriteArray(b7_ds.GetRasterBand(1).ReadAsArray())
# 5. удаляя созданную переменную мы заставляем GDAL записать датасет в файл Composit_rgb.TIF на диск
# Чистим оперативную память
del new_dataset
# ПРОСМОТР ЧЕРЕЗ ArcGis или QGIS
Извлекаем из полученного растра видимый красный и IR и рассчитываем индекс NDVI.
landsat_ds = gdal.Open(workdir + r".\Lansat-8_2021_05_16\Composit_bands_1_7.tif")
# Извлешем из многоканального растра видимый красный и ближний ИК каналы
red_band_array = landsat_ds.GetRasterBand(4).ReadAsArray()
nir_band_array = landsat_ds.GetRasterBand(5).ReadAsArray()# Расчет индекса растительности NDVI
ndvi_array = (nir_band_array - red_band_array) / (nir_band_array + red_band_array)
# Нормируем к единице
ndvi_array_norm = ndvi_array / numpy.nanmax(ndvi_array)
plt.close()
plt.imshow(ndvi_array_norm)
plt.colorbar()
plt.show()
Вот еще несколько интересных индексов, которые также рассматриваются для решения наших задач. Далее расчеты производить не буду, я думаю их механизм ясен и он однотипен, просто приведу расчетные формулы и интерпретации индексов.
Стандартизованный индекс различий застройки (NDBI) | Использует каналы NIR и SWIR для выделения областей застройки. Этот коэффициент позволяет приглушать разницу в освещении поверхности, а также атмосферные эффекты. | NDBI = (SWIR - NIR) / (SWIR + NIR), где SWIR-значения пикселов из коротковолнового инфракрасного канала, NIR-значения пикселов из ближнего инфракрасного канала |
Индекс выгоревших областей (BAI) | Использует значения отражения в красной и ближней инфракрасной области спектра для идентификации областей поверхности, подвергшихся огню. | BAI = 1/((0.1 -RED)^2 + (0.06 - NIR)^2), где Red-значения пикселов из красного канала, NIR-значения пикселов из ближнего инфракрасного канала |
Модифицированный стандартизованный индекс различий воды (MNDWI) | Используется для улучшения отображения объектов открытых водных пространств. Он также снижает значения областей застройки, которые часто коррелированы с открытыми водными пространствами в других индексах. | MNDWI = (Green - SWIR) / (Green + SWIR), где Green-значения пикселов из зеленого канала, SWIR-значения пикселов из коротковолнового инфракрасного канала |
Железистые минералы. Коэффициент железистых минералов (Ferrous Minerals Ratio) | Выделяет все железосодержащие материалы. Он использует соотношение между каналом SWIR и каналом NIR. | Ferrous Minerals Ratio = SWIR / NIR, где SWIR-значения пикселов из коротковолнового инфракрасного канала, NIR-значения пикселов из ближнего инфракрасного канала |
Стандартизованный индекс различий снежного покрова (NDSI) | Разработан для использования данных MODIS (каналы 4 и 6) и Landsat TM (каналы 2 и 5) с целью идентификации снежного покрова при игнорировании облачного покрова. | NDSI = (Green - SWIR) / (Green + SWIR), где Green-значения пикселов из зеленого канала, SWIR-значения пикселов из коротковолнового инфракрасного канала |
Таки м образом, логика веб-ГИС-приложения уже реализует отдельные элементы обработки и анализа полученных растров. Идем дальше, будет интересно, подписывайтесь.
