Эта статья является продолжением цикла статей посвященных развитию стартапа "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()

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

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