Космическая съёмка Земли


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

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

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

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

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


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


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

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


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

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

    Характеристики съёмочный аппаратуры WorldView-3
    Спектральный диапазон Наименование Диапазон, нм Пространственное разрешение в надире, м
    Панхроматический (1 канал; охватывает видимую часть спектра) 450-800 0,31
    Мультиспектральный (8 каналов) Прибрежный 400-450 1,24
    Синий 450-510
    Зелёный 510-580
    Желтый 585-625
    Красный 630-690
    Красный край 705-745
    Ближний инфракрасный 770-895
    Ближний инфракрасный 860-1040
    Многополосный в коротком инфракрасном диапазоне (8 каналов) SWIR-1 1195-1225 3,70
    SWIR-2 1550-1590
    SWIR-3 1640-1680
    SWIR-4 1710-1750
    SWIR-5 2145-2185
    SWIR-6 2185-2225
    SWIR-7 2235-2285
    SWIR-8 2295-2365

    Помимо приведённых каналов 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татистически.

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

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

    $B_\lambda=K_\lambda\;\ast\;DN_\lambda+\;C_\lambda,$


    где $B_\lambda$ — энергетическая яркость для спектральной зоны $\lambda$;
    $DN_\lambda$ — сырые значения яркости;
    $K_\lambda$ — калибровочный коэффициент;
    $C_\lambda$ — калибровочная константа.

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


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

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


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

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


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

    Примеры индексных изображений
    Название индекса Формула Применение
    Индекс содержания оксида железа Красный/синий Для выявления содержания окислов железа
    Индекс содержания глинистых минералов Отношение значений яркости в пределах среднего инфракрасного канала (CИК). CИК1 / CИК2, где CИК1 это диапазон от 1,55 до 1,75 мкм, CИК2 это диапазон от 2,08 до 2,35 мкм Для выявления содержания глинистых минералов
    Индекс содержания железистых минералов Отношение значения яркости в среднем инфракрасном (СИК1; от 1,55 до 1,75 мкм) канале к значению яркости в ближнем инфракрасном канале (БИК). СИК1 / БИК Для выявления содержания железистых минералов
    Индекс красноцветности (RI) Основан на различии отражательной способности красноцветных минералов в красном (К) и зеленом (З) диапазонах. RI = (К — З) / (К + З) Для выявления содержания оксида железа в почве
    Нормализованный дифференциальный индекс снега(NDSI) NDSI – это относительная величина, характеризуемая различием отражательной способности снега в красном (К) и коротковолновом инфракрасном (КИК) диапазонах.
    NDSI=(К — КИК) / (К + КИК)
    Используется для выделения территорий, покрытых снегом. Для снега NDSI > 0,4
    Водный индекс (WI) WI=0,90 мкм / 0,97 мкм Применяется для определения содержания воды в растительности по гиперспектральным снимкам
    Нормализованный дифференциальный вегетационный индекс (NDVI) Хлорофилл листьев растений отражает излучение в ближнем инфракрасном (БИК) диапазоне электромагнитного спектра и поглощает в красном (К). Отношение значений яркости в этих двух каналах позволяет четко отделять и анализировать растительные от прочих природных объектов.
    NDVI=(БИК — К) / (БИК + К)
    Показывает наличие и состояние растительности. Значения NDVI варьируют в пределах от -1 до 1;
    Густая растительность: 0,7;
    Разряженная растительность: 0,5;
    Открытая почва: 0,025;
    Облака: 0;
    Снег и лед: -0,05;
    Вода: -0,25;
    Искусственные материалы (бетон, асфальт): -0,5

    Работа со спутниковыми снимками на 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
    
    #читаем снимок и geojson файл с полигонами
    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.

    image
    Max-pooling

    image
    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.

    Список источников

    Tinkoff.ru 153,26
    Самый большой онлайн-банк в мире
    Поделиться публикацией

    Вакансии компании Tinkoff.ru

    Комментарии 25
    • +1
      Здравствуйте. Я недавно пытался разобраться с кластеризацией объектов на мультиспектральном снимке используя и цвет, и текстуру; не нашёл ни одной метрики, которая бы учитывала одновременно и то, и другое. Не подскажете, это я плохо искал, или действительно задачу ещё не решили?
      • 0
        не понял до конца вопроса. попросил уточнить в лс
      • +1
        Неплохой материал и хорошо что указали источники.
        • +1
          Кто тут говорит, что комментария бесполезны и переменные говорят сами за себя. Но вот авторы пошли дальше, тут даже комментария врут.
          #читаем shape файл с полигонами
          shapefile = fiona.open('buildings.geojson', "r")
          geoms = [feature["geometry"] for feature in shapefile]
          

          • 0
            совсем не понял на счёт комментария
            • 0
              Разница между ERSI Shape форматом и geojson огромна.
              Т.е. название переменной говорит об одном, а по факту там совсем другое. Это как если бы там было написано:
              # открываем вордовский документ
              microsoftWordDocument = open('readme.txt', 'r')
              text = microsoftWordDocument.read()
              
              • 0
                Ага, поправил. Спасибо.
          • +2
            Всё же мучает меня вопрос — зачем кредитной организации: «необходимо разработать алгоритм, выполняющий сегментацию зданий, автомобилей и считающий количество автомобилей на спутниковых снимках»? Это для актуализации данных по должникам?
            И ещё: конкурс на создание алгоритма объявили чтобы не оплачивать его разработку, а получить «на халяву» в рамках мероприятия?
            • 0
              В начале статьи написано для чего — в продуктах для предпринимательства и инвестиций.

              Проведя соревнование, мы не получим «халявное решение». Собственные наработки уже есть. А данных, что подготовлены для соревнования, недостаточно для разработки алгоритма, который бы мог использоваться в реальном продукте.
              • 0
                В свое время страховые проверяли как засевали площади обратившиеся к ним фермеры. Оказались даже не засевали — а сразу осеньюу бежали просить компенсации за убытки. Когда сравни айрофотосьемки весной и летом — оказались поля были не тронуты.
                • 0
                  Про то же подумал. Ещё про строительство и покупку техники. Удивительно, что этим должны заниматься банки и страховые, почему нельзя заказать обработку у профессионалов, например, в Институте вычислительной техники в Новосибе, или в каком-нибудь из осколков ВАСХНИЛ, у экологов из МГУ, у гидрогеологов, в военном вузе… тема анализа мультиспектральных изображений много где живет. Видимо, местный Сам решил, что ИТ-отдел убыточен и вообще бездельники, и надо идти в ногу со времен — нанотехнологии, скрам, блокчейн.
                  • 0
                    Только не у экологов в МГУ — там все ужасно довольно. Про остальные — не знаю.
                    Владение собственной реализацией общей технологии вполне может иметь смысл. Другое дело, что обоснования конкретной формулировки мы, вероятнее всего, не услышим.
                    • 0
                      которые с кафедры экологии — там да, нерадужно, а вот люди, которые разрабатывали дефолиант и антидефолиант по мотивам Вьетнама, до сих пор продуктивно работают на высших и низших травках.
                      А по поводу источников — вбейте в Киберленинке «мультиспектральное изображение» и гляньте афиляции. Неожиданно подборка достаточно свежая, полная и понятная, я оттуда почерпнул больше, чем из англоязычных источников.
              • +1
                Также стоит упомянуть цветовую коррекцию спутниковых снимков — приведение снимка в натуральных цветах (красный, зелёный, синий) к более привычному для человека виду.

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

                Что же касается определения высоты зданий по снимкам высокого разрешения (не стерео) — это, в общем случае, довольно дохлый номер, так как задача имеет сколько-нибудь точное решение только в ситуации, когда здание стоит на ровной горизонтальной поверхности, на которой не растут деревья. В любом другом случае, решение будет гаданием на кофейной гуще. Тем более непонятно, зачем нужна конкретно эта задача, если существуют данные радарной съемки высокого разрешения, обработка которых куда менее однозначна.
                • 0
                  Читать «куда менее неоднозначна» в последнем абзаце.
                  • 0
                    Идея считать высоты зданий по теням слишком красивая, чтобы так просто от неё отказаться. Даже с учётом того, что мелкие здания не получится найти и наличием погрешности для высоких. Думаю стоит по крайней мере попытаться реализовать идею.
                    • 0
                      Красивая — это когда данные для этого подходят.
                • +4
                  Вот эта схема
                  image
                  содержит фундаментальную ошибку: горизонтальная сторона треугольника sl перпендикулярна плоскости стены здания, тогда как видно, что боковая сторона тени лежит к этой плоскости под углом. Таким образом, сторона sl должна быть, на самом деле, длиннее. Это все неплохо вписывается в ситуацию с высокой степенью неточности подобных измерений.
                  • 0
                    Модель сильно упрощённая. Согласен есть недостатки. Можно ещё, например, учитывать satellite elevation angle.
                    • 0
                      Угол будет важен только при съёмке с существенным отклонением камеры от линии надира. Тогда можно замерить видимую высоту здания в пикселах и далее действовать по теореме Пифагора. При условии, что у здания вертикальные стены. Это — апостериорное знание.

                      При горизонтальной почве (парковки, газоны) вокруг здания необходимо отличать сму тень от тёмной стены, примыкающей к ней (снимки с вертикальной камера практически не встречаются) для вычисления высоты кромки крыши.
                  • +3
                    В ГСК единицы измерения угловые и координаты представлены десятичными градусами. Наиболее известная ГСК — WGS 84 (другое название EPSG:4326).

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

                    Возможность линейной интерполяции — не главное свойство систем координат, спроецированных на плоскость. Географические координаты, к слову, тоже можно линейно интерполировать, только вопрос — что из этого получится? (Это касается также и множества плоских систем.) А фраза о том, что наиболее известная (кому, согласно какой статистике?) система координат — это проекция Меркатора — вообще за гранью, потому что проекция — это часть системы координат, способ ее отображения на плоскоть, она не эквивалентна всей системе.

                    Что-то как-то попахивает безграмотностью и верхоглядством.
                    • 0
                      Если сделаем линейную интерполяцию для географических координат, то результат получиться не очень:) Имелось ввиду, что зная геокоординаты объекта и перейдя к ПСК, можно найти координаты объекта уже на самом изображении.

                      Вот здесь спасибо, что заметили неточность. Поправлю. Как-то смешал проекцию и результат применения проекции. Проекция Меркатора — одна из наиболее известных проекций. Хотя бы потому что разработана самим Меркатором, который жил в 16 веке.

                      И комментарий ниже на счёт 37-й зоны справедливый. Всё не раскрыл в статье.

                      Оставлю ссылку на ваши статьи: ликбез по картографическим проекциям habrahabr.ru/post/235283 и Google Web Mercator habrahabr.ru/post/239251. В них довольно подробно всё описано.
                      • 0
                        Меж тем, куча задач вполне себе решаются в географической проекции, то есть ровно в ситуации интерпретации географических (угловых) координат, как прямоугольных декартовых. Результат, и правда, «не очень», но это, нередко, все же делается.

                        Известность — субъективный критерий. Тем более, что если спросить тех, кто эту проекцию использует, не все даже знают, что Меркатор — это фамилия (точнее — латинизированный псевдоним), а не просто название. Используемость — куда более объективный. Вот тут — да, проекция Меркатора — одна из самых часто используемых в веб-приложениях. Это вопрос, скорее, к стилистике изложения и к тому, чтобы не использовать клише вроде «самый известный».
                    • +3
                      Еще, в примерах используется epsg:32637, и ни слова о том, почему выбрана проекция UTM (а не «самая известная» проекция Меркатора, например) и почему — именно 37-й зоны. А то ведь кто-нибудь из Екатеринбурга или Владивостока возьмет, и на свою территорию тоже так сделает :)
                      • 0
                        Спасибо за статью, особенно за таблицу с примерами индексов.
                        На любительском уровне интересовался данной темой, пытаясь найти способ определения качества/возраста асфальтового покрытия. В идеале хотел идентифицировать грунтовое покрытие, каменистое (гравий) и свежеуложенный/старый асфальт. К сожалению, безуспешно.
                        Нашел одну работу по этой теме Monitoring Asphalt Pavement Damages using Remote Sensing Techniques, но моих (не)знаний не хватило чтобы сделать нечто рабочее. Возможно из-за недостаточного разрешения Sentinel-2 в большинстве диапазонов.

                        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                        Самое читаемое