Космическая съёмка Земли23.11.2017 14:03

Cпутниковый снимок в ложных цветах (зелёный, красный, ближний инфракрасный) с пространственным разрешением 3 метра и наложенной маской зданий из OpenStreetMap (спутниковая группировка PlanetScope)

Привет, Хабр! Мы постоянно расширяем источники данных, которые используем для аналитики, поэтому решили добавить ещё и спутниковые снимки. У нас аналитика по спутниковым снимкам полезна в продуктах для предпринимательства и инвестиций. В первом случае статистика по геоданным поможет понять, в каком месте стоит открывать торговые точки, во втором позволяет анализировать деятельность компаний. Например, для строительных компаний можно посчитать, сколько за месяц было построено этажей, для сельскохозяйственных компаний — сколько гектаров урожая взошло и т.д.

В этой статье я постараюсь дать примерное представление о космической съёмке Земли, расскажу о трудностях, с которыми можно столкнуться, начиная работу со спутниковыми снимками: предварительная обработка, алгоритмы для анализа и библиотеки Python для работы со спутниковыми снимками и геоданными. Так что все, кому интересна область компьютерного зрения, добро пожаловать под кат!

Итак, начнём с областей спектра, в которых спутники выполняют съёмку, и рассмотрим характеристики съёмочной аппаратуры некоторых спутников.

Спектральные сигнатуры, атмосферные окна и спутники

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

Отражательная способность сухой травы, асфальта и свежой базальтовой лавы

Как и поверхность Земли, газы в атмосфере также имеют уникальные спектральные сигнатуры. При этом не каждое излучение проходит через атмосферу Земли. Спектральные диапазоны, проходящие через атмосферу, называют «атмосферными окнами», а спутниковые сенсоры настраиваются для измерения в этих окнах.

Атмосферные окна

Рассмотрим спутник с одним из самых высоких пространственных разрешений — WorldView-3 (Пространственное разрешение — величина, характеризующая размер наименьших объектов, различимых на изображении; надир — направление, совпадающее с направлением действия силы гравитации в данной точке).

Помимо приведённых каналов WorldView-3 имеет ещё 12 каналов, предназначенных специально для атмосферной коррекции, — CAVIS (Clouds, Aerosols, Vapors, Ice, and Snow) с разрешением 30 м в надире и длинами волн от 0,4 до 2,2 мкм.

Пример снимка с WorldView-2; панхроматический канал

Пример снимков с различным пространственным разрешением

Другие интересные спутники — SkySat-1 и его близнец SkySat-2. Интересны они тем, что умеют снимать видео продолжительностью до 90 секунд  над одной территорией с частотой 30 кадров/сек и пространственным разрешением 1,1 м.

Видеосъёмка со спутников SkySat

Пространственное разрешение панхроматического канала — 0,9 м, мультиспектральных каналов (синий, зелёный, красный, ближний инфракрасный) — 2 м.

Ещё некоторые примеры:

1. Группировка спутников PlanetScope выполняет съёмку с пространственным разрешением 3 м в красном, зелёном, синем и ближнем инфракрасном диапазоне;
   
2. Группировка спутников RapidEye выполняет съёмку с пространственным разрешением 5 м в красном, крайнем красном, зелёном, синем и ближнем инфракрасном диапазоне;
   
3. Серия российских спутников «Ресурс-П» выполняет съёмку с пространственным разрешением 0,7–1 м в панхроматическом канале, 3–4 м в мультиспектральных каналах (8 каналов).
   

Гиперспектральные сенсоры в отличие от мультиспектральных делят спектр на множество узких диапазонов (примерно 100–200 каналов) и имеют пространственное разрешение уже другого порядка: 10–80 м.

Гиперспектральные снимки доступны не так широко, как мультиспектральные. Космических аппаратов, на борту которых установлены гиперспектральные сенсоры, немного. Среди них Hyperion на борту спутника NASA EO-1 (выведен из эксплуатации), CHRIS на борту спутника PROBA, принадлежащего Европейскому космическому агентству, FTHSI на борту спутника MightySatII исследовательской лаборатории военно-воздушных сил США, ГСА (гиперспектральная аппаратура) на российских космических аппаратах «Ресурс-П».

Обработанный гиперспектральный снимок с бортового сенсора CASI 1500; до 228 каналов; спектральный диапазон 0,4 — 1 нм

Обработанный снимок с космического сенсора EO-1; 220 каналов; спектральный диапазон 0,4 — 2,5 нм

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

Предварительная обработка снимков

Перед тем, как приступить к анализу снимков, необходимо выполнить их предварительную обработку:

1. Геопривязка;
   
2. Ортокоррекция;
   
3. Радиометрическая коррекция и калибровка;
   
4. Атмосферная коррекция.
   

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

Также стоит упомянуть цветовую коррекцию спутниковых снимков — приведение снимка в натуральных цветах (красный, зелёный, синий) к более привычному для человека виду.

Приблизительная геопривязка вычисляется по исходному положению спутника на орбите и геометрии изображения. Уточнение геопривязки выполняется по наземным точкам привязки (Ground Control Points — GCP). Эти контрольные точки ищутся на карте и на снимке, а зная их координаты в разных системах координат, можно найти параметры преобразования (конформное, аффинное, перспективное либо полиномиальное) из одной системы координат в другую. Поиск GСP выполняется при помощи GPS-съёмки [ист. 1 с. 230, ист. 2] либо при помощи сопоставления двух снимков, на одном из которых точно известны координаты GCP, по ключевым точкам.

Ортокоррекция изображения — процесс геометрической коррекции изображений, при котором устраняются перспективные искажения, развороты, искажения вызванные дисторсией объектива и другие. Изображение при этом приводится к плановой проекции, то есть такой, при которой каждая точка местности наблюдается строго вертикально, в надир.

Так как спутники осуществляют съемку с очень большой высоты (сотни километров), то при съёмке в надире искажения должны быть минимальными. Но космический аппарат не может всё время снимать в надире, иначе пришлось бы очень долго ждать момента,  когда он пройдет над заданной точкой. Для устранения этого недостатка спутник «доворачивают», и большинство кадров получаются перспективными. Следует заметить, что углы съемки могут достигать 45 градусов, и при большой высоте это приводит к значительным искажениям.

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

Модель камеры спутника представляется в виде обобщенных аппроксимирующих функций (рациональных полиномов — RPC коэффициентов), а высотные данные могут быть получены в результате наземных измерений, при помощи горизонталей с топографической карты, стереосъемки, по радарным данным или из общедоступных грубых цифровых моделей рельефа: SRTM (разрешение 30–90 м) и ASTER GDEM (разрешение (15–90 м).

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

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

Удаление сбойных пикселей и вертикальных полос

Радиометрическая коррекция выполняется двумя методами:

1. C использованием известных параметров и настроек съемочного прибора (корректировочных таблиц);
   
2. Cтатистически.
   

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

Радиометрическая калибровка снимков — перевод «сырых значений» яркости в физические единицы, которые можно сопоставлять с данными других снимков:

где  — энергетическая яркость для спектральной зоны ;
 — сырые значения яркости;
 — калибровочный коэффициент;
 — калибровочная константа.

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

Факторы, влияющие на попадание отраженной солнечной радиации на сенсоры спутника

Существуют различные алгоритмы выполнения атмосферной коррекции (например метод DOS — Dark Object Subtraction). Входными параметрами для моделей служат: геометрия расположения Солнца и датчика, атмосферная модель для газообразных компонентов, модель аэрозоля (тип и концентрация), оптическая толщина атмосферы, коэффициент поверхностного отражения и спектральные каналы.
Для атмосферной коррекции можно также применять алгоритм удаления дымки с изображения — Single Image Haze Removal Using Dark Channel Prior (реализация).

Пример работы Single Image Haze Removal Using Dark Channel Prior

Индексные изображения

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

Работа со спутниковыми снимками на Python

Одним из форматов, в которых принято хранить спутниковые снимки, является GeoTiff (ограничимся только им). Для работы с GeoTiff на Python можно воспользоваться библиотекой gdal либо rasterio.

Для установки gdal и rasterio лучше воспользоваться conda:

conda install -c conda-forge gdal
conda install -c conda-forge rasterio

Другие библиотеки для работы со спутниковыми снимками легко ставятся через pip.

Чтение GeoTiff через gdal:

from osgeo import gdal
import numpy as np

src_ds = gdal.Open('image.tif')
img = src_ds.ReadAsArray()
#height, width, band
img = np.rollaxis(img, 0, 3)

width = src_ds.RasterXSize
height = src_ds.RasterYSize

gt = src_ds.GetGeoTransform()
#геокоординаты крайних точек снимка
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 
maxx = gt[0] + width*gt[1] + height*gt[2]
maxy = gt[3]

Систем координат у спутниковых снимков довольно много. Их можно разделить на две группы: географические системы координат (ГСК) и плоские системы координат (ПСК).

В ГСК единицы измерения угловые и координаты представлены десятичными градусами. Наиболее известная ГСК — WGS 84 (другое название EPSG:4326).

В ПСК единицы измерения линейные и координаты могут быть выражены как метры, футы, километры и т.д., поэтому их можно линейно интерполировать. Наиболее известная ПСК — проекция Меркатора.

Карту (разметку снимка) принято хранить не в растровом виде, а в виде точек, линий и полигонов. Внутри файлов с разметкой снимков хранятся геокоординаты вершин этих геометрических объектов. Для чтения и работы с ними можно использовать библиотеки fiona и shapely.

Скрипт для растеризации полигонов:

import numpy as np
import rasterio
#fiona умеет работать с разными форматами хранения полигонов, например shp и geojson
import fiona
import cv2

from shapely.geometry import shape
from shapely.ops import cascaded_union
from shapely.ops import transform
from shapely.geometry import MultiPolygon

def rasterize_polygons(img_mask, polygons):
    if not polygons:
        return img_mask
    
    if polygons.geom_type == 'Polygon':
        polygons = MultiPolygon([polygons])
    
    int_coords = lambda x: np.array(x).round().astype(np.int32)
    exteriors = [int_coords(poly.exterior.coords) for poly in polygons]
    interiors = [int_coords(pi.coords) for poly in polygons
                 for pi in poly.interiors]
    
    cv2.fillPoly(img_mask, exteriors, 255)
    cv2.fillPoly(img_mask, interiors, 0)
    return img_mask


def get_polygons(shapefile):
    
    geoms = [feature["geometry"] for feature in shapefile]
    polygons = list()
    for g in geoms:
            s = shape(g)
            if s.geom_type == 'Polygon':
                if s.is_valid:
                    polygons.append(s)
                else:
                    #устраняем самопересечения полигона
                    polygons.append(s.buffer(0))
            elif s.geom_type == 'MultiPolygon':
                for p in s:
                    if p.is_valid:
                        polygons.append(p)
                    else:
                        #устраняем самопересечения полигона
                        polygons.append(p.buffer(0))

    mpolygon = cascaded_union(polygons)
    
    return mpolygon

#читаем снимок и shape файл с полигонами
src = rasterio.open('image.tif')
shapefile = fiona.open('buildings.geojson', "r")

#геокоординаты снимка (системы координат снимка и полигонов должны совпадать)
left = src.bounds.left
right = src.bounds.right
bottom = src.bounds.bottom
top = src.bounds.top

height,width = src.shape

#извлекаем мультиполигоны
mpolygon = get_polygons(shapefile)

#так как система координат плоская, то можем получить координаты полигонов уже на самом изображении
mpolygon = transform(lambda x, y, z=None: (width * (x - left) / float(right - left),\
                                                       height - height * (y - bottom) / float(top - bottom)), mpolygon)

#рисуем маску
real_mask = np.zeros((height,width), np.uint8)
real_mask = rasterize_polygons(real_mask, mpolygon)

Во время дешифрования снимков может понадобиться операция пересечения полигонов (например, автоматически разметив здания по карте, хочется также автоматически убрать разметку зданий в тех местах снимка, где есть облака). Для этого есть алгоритм Уайлера-Атертона, но он работает только с полигонами без самопересечений. Для устранения самопересечений нужно проверить пересечение всех ребер с другими ребрами полигона и добавить новые вершины. Эти вершины разобьют соответствующие им ребра на части. В библиотеке shapely есть метод для устранения самопересечений — buffer (0).

Для перевода из ГСК в ПСК можно воспользоваться библиотекой PyProj (либо сделать это в rasterio):

#переводим координаты из географической системы координат epsg:4326 в координаты плоской системы координат epsg:32637 
def wgs84_to_32637(lon, lat):
    from pyproj import Proj, transform
    
    inProj = Proj(init='epsg:4326')
    outProj = Proj(init='epsg:32637')
    
    x,y = transform(inProj,outProj,lon,lat)
    
    return x,y

Метод главных компонент

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

Такое преобразование также проводят для серии разновременных снимков, приведенных в единую систему координат, для выявления динамики, которая ярко проявляется в одной или двух компонентах.

Скрипт для сжатия 4-х канального изображения до 3-х канального:

from osgeo import gdal
from sklearn.decomposition import IncrementalPCA
from sklearn import preprocessing
import numpy as np
import matplotlib.pyplot as plt

#читаем снимок, меняем порядок на height, width, band и вырезаем часть изображения, чтобы PCA быстрее справился с данными
src_ds = gdal.Open('0.tif')
img = src_ds.ReadAsArray()
img = np.rollaxis(img, 0, 3)[300:1000, 700:1700, :]
height, width, bands = img.shape

#готовим данные для PCA
data = img.reshape((height * width, 4)).astype(np.float)
num_data, dim = data.shape
pca = IncrementalPCA(n_components = 3, whiten = True)

#cтандартизация данных (достаточно и центрирования)
x = preprocessing.scale(data)

#сжимаем каналы и приводим к виду удобному для визуализации
y = pca.fit_transform(x)
y = y.reshape((height, width, 3))

#нормализуем яркость пикселей
for idx in range(0, 3):
    band_max = y[:, :, idx].max()
    y[:, :, idx] = np.around(y[:, :, idx] * 255.0 / band_max).astype(np.uint8)

Снимок с группы спутников PlanetScope (red, green, blue без цветовой коррекции)

Снимок с группы спутников PlanetScope (green, red, near infrared)

Снимок, полученный при помощи метода главных компонент

Метод спектрального разделения (Spectral Unmixing)

Метод спектрального разделения применяют для распознавания на снимках объектов, размер которых значительно меньше размера пикселя.

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

Смешанная спектральная кривая

Дешифрованный снимок

Сегментация спутниковых снимков

На данный момент state-of-the-art результаты в задачах бинарной сегментации изображений показывают модификации U-Net модели.

Архитектура U-Net модели (размер изображения на выходе меньше размера входного изображения; делается это из-за того, что сеть хуже предсказывает на краях изображения)

Автор U-Net разработал архитектуру на основе другой модели — Fully convolutional network (FCN), особенностью которой является наличие только свёрточных слов (не считая max-pooling).
U-Net отличается от FCN тем, что добавлены слои, в которых max-pooling заменён на up-convolution. Таким образом, новые слои постепенно увеличивают выходное разрешение. Также признаки с энкодер-части комбинируются с признаками из декодер-части, чтобы модель могла делать более точные предсказания за счёт дополнительной информации.

Модель, в которой отсутствует проброс признаков из энкодер-части в декодер-часть, называется SegNet и на практике показывает результаты хуже U-Net.

Max-pooling

Up-convolution

Отсутствие в U-Net, Segnet и FCN слоёв, привязанных к размеру изображения, позволяет подавать на вход одной и той же сети изображения разного размера (размер снимка должен быть кратным количеству фильтров в первом свёрточном слое).

В keras это реализуется так:

inputs = Input((channel_number, None, None))

Обучение, как и предсказание, можно проводить либо на фрагментах изображения (кропах), либо на всём изображении целиком, если позволяет память GPU. При этом в первом случае:
1) больше размер батчей, что хорошо скажется на точности модели, если данные зашумлены и неоднородны;
2) меньше риск переобучения, т.к. данных гораздо больше, чем при обучении на изображениях полного размера.

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

U-Net  —  простая и мощная архитектура для задачи бинарной сегментации, на github можно найти не одну реализацию для любого DL фреймворка, но для сегментации большого количества классов эта архитектура проигрывает другим архитектурам, например, PSP-Net. Здесь можно прочитать интересный обзор по архитектурам для семантической сегментации изображений.

Определение высоты зданий

Высоту зданий можно определять по их теням. Для этого необходимо знать: размер пикселя в метрах, длину тени в пикселях и sun (solar) elevation angle (угол солнца над горизонтом).

Геометрия солнца, спутника и здания

Вся сложность задачи в том, чтобы как можно точнее сегментировать тень здания и определить длину тени в пикселях. Проблем добавляет ещё и наличие облаков на снимках.

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

import pandas as pd
import numpy as np
import rasterio
import math
import cv2

from shapely.geometry import Point

#фильтр для бинарных изображений по размеру сегментов (реализован на основе волнового алгоритма)
#лучше воспользоваться функциями remove_small_objects и remove_small_holes из библиотеки skimage
def filterByLength(input, length, more=True):
    import Queue as queue
    
    copy = input.copy()
    output = np.zeros_like(input)
    
    for i in range(input.shape[0]):
        
        for j in range(input.shape[1]):
            if (copy[i][j] == 255):
                 
                q_coords = queue.Queue()
                output_coords = list()
                
                copy[i][j] = 100
                
                q_coords.put([i,j])
                output_coords.append([i,j])
                
                while(q_coords.empty() == False):
                    
                    currentCenter = q_coords.get()
                    
                    for idx1 in range(3):
                        for idx2 in range(3):
                            
                            offset1 = - 1 + idx1
                            offset2 = - 1 + idx2
                            
                            currentPoint = [currentCenter[0] + offset1, currentCenter[1] + offset2]

                            if (currentPoint[0] >= 0 and currentPoint[0] < input.shape[0]):
                                    if (currentPoint[1] >= 0 and currentPoint[1] < input.shape[1]):
                                        if (copy[currentPoint[0]][currentPoint[1]] == 255):

                                            copy[currentPoint[0]][currentPoint[1]] = 100

                                            q_coords.put(currentPoint)
                                            output_coords.append(currentPoint)
                
                if (more == True):
                    if (len(output_coords) >= length):
                        for coord in output_coords:
                            output[coord[0]][coord[1]] = 255
                else:
                    if (len(output_coords) < length):
                        for coord in output_coords:
                            output[coord[0]][coord[1]] = 255
                
    return output

#переводим координаты из плоской системы координат epsg:32637 в координаты географической системы координат epsg:4326
def getLongLat(x1, y1):
    from pyproj import Proj, transform
    
    inProj = Proj(init='epsg:32637')
    outProj = Proj(init='epsg:4326')
    
    x2,y2 = transform(inProj,outProj,x1,y1)
    
    return x2,y2

#переводим координаты из географической системы координат epsg:4326 в координаты плоской системы координат epsg:32637 
def get32637(x1, y1):
    from pyproj import Proj, transform
    
    inProj = Proj(init='epsg:4326')
    outProj = Proj(init='epsg:32637')
    
    x2,y2 = transform(inProj,outProj,x1,y1)
    
    return x2,y2

#переводим координаты из плоской системы координат epsg:32637 в координаты изображения(пиксели)
def crsToPixs(width, height, left, right, bottom, top, coords):
    
    x = coords.xy[0][0]
    y = coords.xy[1][0]
    
    x = width * (x - left) / (right - left)
    y = height - height * (y - bottom) / (top - bottom)
    
    x = int(math.ceil(x))
    y = int(math.ceil(y))
    
    return x,y

#сегментация тени при помощи порогового преобразования
def shadowSegmentation(roi, threshold = 60):
    
    thresh = cv2.equalizeHist(roi)
    ret, thresh = cv2.threshold(thresh,threshold,255,cv2.THRESH_BINARY_INV)
    
    tmp = filterByLength(thresh, 50)
    
    if np.count_nonzero(tmp) != 0:
        thresh = tmp
        
    return thresh

#определяем размер(длину) тени; x,y - координаты здания на изображении thresh
def getShadowSize(thresh, x, y):
    
    #определяем минимальную дистанцию от здания до пикселей тени 
    min_dist = thresh.shape[0]
    min_dist_coords = (0, 0)
    
    for i in range(thresh.shape[0]):
        for j in range(thresh.shape[1]):
            if (thresh[i,j] == 255) and (math.sqrt( (i - y) * (i - y) + (j - x) * (j - x) ) < min_dist):
                min_dist = math.sqrt( (i - y) * (i - y) + (j - x) * (j - x) )
                min_dist_coords = (i, j) #y,x
                
    #определяем сегмент, который содержит пиксель с минимальным расстоянием до здания
    import Queue as queue
    
    q_coords = queue.Queue()
    q_coords.put(min_dist_coords)
    
    mask = thresh.copy()
    output_coords = list()
    output_coords.append(min_dist_coords)
    
    while q_coords.empty() == False:
        currentCenter = q_coords.get()
        
        for idx1 in range(3):
            for idx2 in range(3):

                offset1 = - 1 + idx1
                offset2 = - 1 + idx2

                currentPoint = [currentCenter[0] + offset1, currentCenter[1] + offset2]

                if (currentPoint[0] >= 0 and currentPoint[0] < mask.shape[0]):
                        if (currentPoint[1] >= 0 and currentPoint[1] < mask.shape[1]):
                            if (mask[currentPoint[0]][currentPoint[1]] == 255):

                                mask[currentPoint[0]][currentPoint[1]] = 100
                                q_coords.put(currentPoint)
                                output_coords.append(currentPoint)
    
    #отрисовываем ближайшую тень
    mask = np.zeros_like(mask)
    for i in range(len(output_coords)):
        mask[output_coords[i][0]][output_coords[i][1]] = 255
    
    #определяем размер(длину) тени при помощи морфологической операции erode
    kernel = np.ones((3,3),np.uint8)
    i = 0
    while np.count_nonzero(mask) != 0:
        mask = cv2.erode(mask,kernel,iterations = 1)
        i += 1
    
    return i + 1

#определяем область, где нет облаков
def getNoCloudArea(b, g, r, n):
    
    gray = (b + g + r + n) / 4.0
    
    band_max = np.max(gray)
    
    gray = np.around(gray * 255.0 / band_max).astype(np.uint8)
    gray[gray == 0] = 255

    ret, no_cloud_area = cv2.threshold(gray, 100, 255, cv2.THRESH_BINARY)
    kernel = np.ones((50, 50), np.uint8)
    no_cloud_area = cv2.morphologyEx(no_cloud_area, cv2.MORPH_OPEN, kernel)
    kernel = np.ones((100, 100), np.uint8)

    no_cloud_area = cv2.morphologyEx(no_cloud_area, cv2.MORPH_DILATE, kernel)
    no_cloud_area = cv2.morphologyEx(no_cloud_area, cv2.MORPH_DILATE, kernel)
    no_cloud_area = 255 - no_cloud_area
    
    return no_cloud_area

#читаем csv в формате lat,long,height
df = pd.read_csv('buildings.csv')
#читаем csv файл с информацией о снимке
image_df = pd.read_csv('geotiff.csv')

#угол солнца над горизонтом
sun_elevation = image_df['sun_elevation'].values[0]

with rasterio.open('image.tif') as src:
    
        #геокоординаты снимка
        #в этом случаем снимок сохранён в ПСК epsg:32637
        left = src.bounds.left
        right = src.bounds.right
        bottom = src.bounds.bottom
        top = src.bounds.top
        
        height,width = src.shape
        
        b, g, r, n = map(src.read, (1, 2, 3, 4))
        
        #берём зелёный канал, т.к. он самый контрастный
        band_max = g.max()
        img = np.around(g * 255.0 / band_max).astype(np.uint8)

        #определяем маску области, где нет облаков
        no_cloud_area = getNoCloudArea(b, g, r, n)

        heights = list()
        
        #определяем тень в окне размером (size, size)
        size = 30

        for idx in range(0, df.shape[0]):

            #геокоординаты и высота здания
            lat = df.loc[idx]['lat']
            lon = df.loc[idx]['long']
            build_height = int(df.loc[idx]['height'])

            #переводим геокоординаты в координаты плоской системы координат epsg:32637
            #(в ПСК можно выполнять линейную интерполяцию и таким образом находить координаты зданий уже на самом изображении)
            build_coords = Point(get32637(lon, lat))
            
            #координаты снимка и зданий в одной ПСК, поэтому можно определить координаты здания уже на самом изображении
            x,y = crsToPixs(width, height, left, right, bottom, top, build_coords)

            #ищем тень зданий, если в этом месте нет облаков 
            if no_cloud_area[y][x] == 255:

                #зная азимут солнца, мы знаем в каком направлении должны быть тени зданий
                #поиск тени зданий будем выполнять в этом направлении
                roi = img[y-size:y,x-size:x].copy()
                shadow = shadowSegmentation(roi)

                #(size, size) - координаты здания в roi
                #умножаем длину тени в пикселях на размер пикселя в метрах
                #(в этом случае пространственное разрешение 3 м)
                shadow_length = getShadowSize(shadow, size, size) * 3

                est_height = shadow_length * math.tan(sun_elevation * 3.14 / 180)        
                est_height = int(est_height)
                
                heights.append((est_height, build_height))

        MAPE = 0
        
        for i in range(len(heights)):
            MAPE += ( abs(heights[i][0] - heights[i][1]) / float(heights[i][1]) )
            
        MAPE *= (100 / float(len(heights)))

Пример работы алгоритма определения высоты зданий по тени

На снимках с группы спутников PlanetScope с пространственным разрешением 3 м ошибка определения высоты зданий по MAPE (mean absolute percentage error) составила ~ 30%. Всего было рассмотрено 40 зданий и один снимок. Однако на субметровых снимках исследователи получили ошибку всего 4–5%.

Заключение

Дистанционное зондирование Земли даёт много возможностей для проведения аналитики. Например, такие компании, как Orbital Insight, Spaceknow, Remote Sensing Metrics, OmniEarth и DataKind, на основе съёмки со спутников выполняют мониторинг производства и потребления в США, анализ урбанизации, трафика, растительности, экономики и т.д. При этом сами снимки становятся всё более доступными. Например, снимки со спутников Landsat-8 и Sentinel-2 с пространственным разрешением больше 10 м находятся в открытом доступе и постоянно обновляются.

В России компании Совзонд, СканЭкс, Ракурс, Гео-Альянс и Северная Географическая Компания также проводят геоаналитику по спутниковым снимкам и являются официальными дистрибьюторами компаний-операторов спутников ДЗЗ: ОАО «Российские космические системы» (Россия), DigitalGlobe (США), Planet (США), Airbus Defence and Space (Франция-Германия) и др.

P.S. Вчера мы запустили онлайн-соревнование по спутниковым снимкам, состоящее из двух задач:

1. Сегментация зданий и автомобилей и подсчёт автомобилей;
   
2. Определение высот зданий.
   

Так что все, кто кому интересна область компьютерного зрения и дистанционного зондирования Земли, присоединяйтесь!

Мы планируем провести ещё одно соревнование по снимкам Роскосмоса, Airbus Defence and Space и PlanetScope.