Удаление облаков со спутниковых изображений Landsat с помощью Python

Расчет облачного покрова в интересующей вас области, удаление облаков и закрашивание их с помощью другого спутникового снимка

Наши технологии продвинулись вперед, чтобы покорить просторы космоса. Мы запустили спутники, оснащенные самыми современными датчиками, для наблюдения за нашими изменяющимися планетами. Тем не менее, один противник унизил эти передовые системы — облака.

Если вы хотите сравнить значения спектральных индексов или обучить модели машинного обучения на спутниковых снимках, облака — это проблема. Итак, мы рассмотрим, как использовать облачную маску, чтобы:

Рассчитайте процент облачности в интересующей вас области.
Удалить облака с изображения
Нарисуйтепиксели облаков, используя изображение, снятое в другой день.

Когда дело доходит до самой маски, мы рассмотрим два варианта:

Тот, который содержится в файле Landsat QA (QA_PIXEL.tif).
Альтернативный подход к машинному обучению

Мы объясним код Python, используемый для этого, и вы можете найти полный проект на GitHub. Давайте очистим эти небеса!

Загрузка спутникового снимка

Обсуждаемые нами методы будут работать с любыми растровыми данными. В этой статье мы специально будем работать со спутниковыми снимками Landsat. Вы можете скачать сцену Landsat с помощью портала EarthExplorer. В качестве альтернативы, если вы хотите использовать Python, статья ниже проведет вас через этот процесс:

Загрузка спутниковых изображений Landsat с помощью Python
Упрощение загрузки сцен Landsat с помощью пакета Python landsatxploretowardsdatascience.com

В итоге у вас должна получиться папка, подобная Рисунок 1. Это все файлы, доступные для научного продукта Landsat level 2. Мы будем работать с выделенными файлами. Это 3 полосы видимого света и файл QA_PIXEL.

Эта конкретная сцена была снята над Кейптауном, Южная Африка. Чтобы увидеть это, мы визуализируем полосы видимого света с помощью функции get_RGB. Это принимает имя файла/идентификатор в качестве параметра. Затем он загрузит бэнды (строки 7–9), сложит их (строка 12), масштабирует (строка 13) и обрежет (строка 16).

import tifffile as tiff
import numpy as np

def get_RGB(ID):

    # Load Blue (B2), Green (B3) and Red (B4) bands
    B2 = tiff.imread('./data/{}/{}_SR_B2.TIF'.format(ID, ID))
    B3 = tiff.imread('./data/{}/{}_SR_B3.TIF'.format(ID, ID))
    B4 = tiff.imread('./data/{}/{}_SR_B4.TIF'.format(ID, ID))

    # Stack and scale bands
    RGB = np.dstack((B4, B3, B2))
    RGB = np.clip(RGB*0.0000275-0.2, 0, 1)

    # Clip to enhance contrast
    RGB = np.clip(RGB,0,0.3)/0.3

    return RGB

Мы используем эту функцию, чтобы получить визуализацию RGB для сцены, которую мы загрузили на рисунке 1 (строки 3–4). Вы можете увидеть полученное изображение на Рис. 2. Обратите внимание на все облака! Посмотрим, сможем ли мы справиться с ними с помощью облачной маски Landsat.

import matplotlib.pyplot as plt

ID = 'LC09_L2SP_175083_20230410_20230412_02_T1'
RGB = get_RGB(ID)

# Plot the RGB image
fig, ax = plt.subplots(figsize=(20, 10))
ax.imshow(RGB)

Облачные маски Landsat

Что такое облачная маска?

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

Существуют различные методы, используемые для классификации облачных пикселей. Landsat использует алгоритм под названием CFMask. При этом используется комбинация различных статистических данных, основанных на спектральных индексах. Он имеет точность классификации 89% для пикселей облаков и 96% для теней от облаков [1]. Позже мы обсудим альтернативный подход машинного обучения к созданию облачных масок.

Понимание диапазона обеспечения качества (QA)

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

Например, предположим, что один из пикселей в файле QA имеет значение 22280. Чтобы понять, что означает это целое число, мы сначала преобразуем его в 16-битное двоичное число — 0101011100001000. Затем каждый бит в этом числе сопоставляется с флагом, показанным на рис. 3. Бит 3 имеет значение 1, указывающее, что этот пиксель классифицируется как облако. В дальнейшем мы будем рассматривать типы покрытия для битов с 1 по 4 для нашей облачной маски.

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

def get_mask(val,type='cloud'):
    
    """Get mask for a specific cover type"""

    # convert to binary
    bin_ = '{0:016b}'.format(val)

    # reverse string
    str_bin = str(bin_)[::-1]

    # get bit for cover type
    bits = {'cloud':3,'shadow':4,'dilated_cloud':1,'cirrus':2}
    bit = str_bin[bits[type]]

    if bit == '1':
        return 0 # cover
    else:
        return 1 # no cover

Чтобы применить эту функцию, мы начинаем с загрузки нашей полосы QA (строка 2). Мы применяем get_mask к каждому пикселю полосы для каждого из 4 типов покрытия (строки 6–9). Вместо того, чтобы перебирать каждый пиксель, мы ускоряем процесс, используя np.vectorize(). Каждая маска будет массивом с теми же размерами, что и сцена, где каждый элемент имеет значение 0 или 1.

# QA band
QA = tiff.imread('./data/{}/{}_QA_PIXEL.TIF'.format(ID, ID))
QA = np.array(QA)

# Get masks
cloud_mask = np.vectorize(get_mask)(QA,type='cloud')
shadow_mask = np.vectorize(get_mask)(QA,type='shadow')
dilated_cloud_mask = np.vectorize(get_mask)(QA,type='dilated_cloud')
cirrus_mask = np.vectorize(get_mask)(QA,type='cirrus')

Мы можем увидеть эти маски в действии на Рис. 4. Здесь мы наложили изображение, которое мы видели на Рисунок 2, с разными цветами в зависимости от классификации в 4 масках. Увеличив масштаб одной области, мы можем лучше понять каждый тип покрытия:

Тени – это более темные области под облаками.
Расширенные облака – это края вокруг пикселей, классифицируемые как облака.
Перистые – "тонкие" облака.

Если пиксель классифицируется как перистый, он также будет классифицироваться как облако.

Приведенный ниже код используется для создания рис. 4. Мы добавляем слой для каждой из масок к исходному изображению RGB (строки 15–22). Это делается с помощью функции cv2.addWeighted() (строка 22). Каждому слою маски присваивается свой цвет (строки 8–11) с прозрачностью 50%. Область интереса (AOI) обрезается путем взятия пикселей в заданном диапазоне (строка 43). Далее мы более подробно рассмотрим эту область.

import cv2
import matplotlib as mpl

# segmentation image
seg = RGB.copy()

# color for each cover type
colors = np.array([[247, 2, 7],
                    [201, 116, 247],
                    [0, 234, 255],
                    [3, 252, 53]])/255

masks = [cloud_mask, shadow_mask, dilated_cloud_mask, cirrus_mask]

for i,mask in enumerate(masks):

    # color for cover type
    temp = seg.copy()
    temp[mask == 0] = colors[i]

    # add to segmentation
    seg = cv2.addWeighted(seg, 0.5, temp, 0.5, 0)

fig, ax = plt.subplots(1,2,figsize=(20, 10))

ax[0].imshow(seg)

# add legend with colors for each cover type
legend_elements = [mpl.patches.Patch(facecolor=colors[0], label='Cloud'),
                     mpl.patches.Patch(facecolor=colors[1], label='Shadow'),
                        mpl.patches.Patch(facecolor=colors[2], label='Dilated Cloud'),
                        mpl.patches.Patch(facecolor=colors[3], label='Cirrus')]

ax[0].legend(handles=legend_elements, loc='upper right')

# draw white rectangle around area of interest
h = 300
x,y = 4500,4500
rect = mpl.patches.Rectangle((x-h,y-h),h*2,h*2,linewidth=2,edgecolor='w',facecolor='none')
ax[0].add_patch(rect)

# crop area of interest
crop_seg =  seg[y-h:y+h,x-h:x+h,:]

ax[1].imshow(crop_seg)

ax[0].set_axis_off()
ax[1].set_axis_off()

Оценка области интереса

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

Чтобы увидеть это, мы выбираем ту же интересующую область, что и на Рисунке 4 (строка 2). Имейте в виду, что облачная маска имеет значение 0 при отсутствии облаков и 1 в противном случае. Мы инвертируем это и берем среднее значение (строка 5). Это дает нам значение 69%. К сожалению, наш AOI более облачен, чем сцена в целом!

# crop area of interest
cloud_aoi = cloud_mask[y-h:y+h,x-h:x+h]

# calculate percentage of cloud cover
per_cloud = np.average(1-cloud_aoi)*100
print('Percentage of clouds: {:.2f}%'.format(per_cloud))

Удаление облаков

Когда у нас есть маска, удалить нежелательные пиксели несложно. Нам просто нужно умножить изображение RGB на маску. Код ниже удалит пиксели, классифицированные как облака (строка 2). Имейте в виду, что изображение RGB представляет собой трехмерный массив. Для сравнения, cloud_mask представляет собой двумерный массив. Вот почему мы используем np.newaxis, чтобы добавить дополнительное измерение.

# Remove clouds
rm_clouds = RGB*cloud_mask[:, :, np.newaxis]

Результат вы можете увидеть на Рис. 5. Все пиксели облаков теперь заменены значениями [0,0,0]. Это то же значение, которое используется для ограничивающей рамки. Это означает, что пиксели не имеют данных и могут быть проигнорированы при дальнейшем анализе. В конечном счете, результаты этого анализа зависят от точности маски облаков. Вот почему исследователи постоянно стремятся улучшить эти маски.

Подход машинного обучения к облачным маскам

Мы упоминали, что Landsat использует алгоритм под названием CFMask. Он основан на дистанционном зондировании и метеорологических знаниях о том, как пиксели облаков представлены на спутниковых изображениях. В качестве альтернативы мы можем использовать машинное обучение. Для этого требуются тренировочные изображения, на которых все пиксели облаков были помечены вручную. Затем алгоритм ML изучает взаимосвязь между спектральными полосами и метками.

Один подход использует генетический алгоритм, обученный на диапазонах прибрежного аэрозоля (B1), красного (B4) и коротковолнового инфракрасного излучения 2 (B7). Результирующий алгоритм задается в функции cloud_pred. Формула pred (строки 4–10) дает взвешивание трех полос. Если взвешенное значение больше полосы B7, пиксель помечен как облако (строка 12).

def cloud_pred(B1,B4,B7):
    """Cloud prediction model"""

    pred = 2.16246741593412 -0.796409165054949*B4 + \
    0.971776520302587*np.sqrt(abs(0.028702220187686*B7*B1 + \
    0.971297779812314*np.sin(B1))) + \
    0.0235599298084993*np.floor(0.995223926146334*np.sqrt(abs(0.028702220187686*B7*B1 + 0.971297779812314*np.sin(B1)))+ \
    0.00477607385366598*abs(0.028702220187686*B7*B1 + \
    0.971297779812314*np.sin(B1))) - 0.180030905136552*np.cos(B4) + \
    0.0046635498889134*abs(0.028702220187686*B7*B1 + 0.971297779812314*np.sin(B1))

    cloud = np.where(pred > B7,0,1)

    return cloud

Чтобы использовать эту функцию, мы начинаем с загрузки 3 необходимых бэндов (строки 4–6). Мы получаем маску облака, передавая эти полосы в функцию cloud_pred (строка 9). Затем мы можем использовать эту маску для удаления облаков так же, как и раньше (строка 12). Полученное изображение показано на Рис. 6.

ID = "LC09_L2SP_175083_20230410_20230412_02_T1"

# Get Coastal Aerosol (B1), Red (B4) and Shortwave Infrared 2 (B7) bands
B1 = tiff.imread('./data/{}/{}_SR_B1.TIF'.format(ID, ID))
B4 = tiff.imread('./data/{}/{}_SR_B4.TIF'.format(ID, ID))
B7 = tiff.imread('./data/{}/{}_SR_B7.TIF'.format(ID, ID))

# Get cloud mask
cloud_mask_2 = cloud_pred(B1,B4,B7)

# Remove clouds
rm_clouds = RGB*cloud_mask_2[:, :, np.newaxis]

Исследователи заявляют, что точность описанного выше подхода составляет 89% — аналогично CFMask. Как уже упоминалось, повышение точности алгоритмов является активной областью исследований. Современные подходы используют алгоритмы глубокого обучения для сегментации изображений, такие как U-Net. Недостатком этих подходов по сравнению с CFMask является необходимость точно размеченных обучающих данных.

Раскрашивание облаков

Для некоторых приложений может быть достаточно просто удалить облака с вашего изображения. Для других мы хотим заменить или «раскрасить» облака. Например, для создания визуализации пейзажа без облаков.

На Рис. 7 мы видим то же облачное изображение, с которым мы работали. Снято 04.10.2023. К счастью, через 8 дней было сделано четкое изображение той же области. Мы рассмотрим, как заменить пиксели в облачном изображении (строка 6) на пиксели в ясном изображении (строка 7).

# IDs for cloudy and non-cloudy images
IDs = ['LC09_L2SP_175083_20230410_20230412_02_T1',
       'LC08_L2SP_175083_20230418_20230429_02_T1']

# Get RGB
cloudy = get_RGB(IDs[0])
clear = get_RGB(IDs[1])

Мы сделаем это как для AOI, так и для всей сцены. В обоих случаях важно учитывать геолокацию пикселей. Несмотря на то, что облачные и четкие изображения сделаны над одной и той же областью, пиксели не будут идеально выровнены. Другими словами, координаты UTM пикселя в позиции массива [x,y] будут разными для каждого изображения. Следовательно, если бы мы полагались только на позиции массива, мы бы в конечном итоге заменили пиксели неправильным участком земли. Чтобы справиться с этим, мы опираемся на методы, используемые в этой статье:

Как наносить координаты на спутниковые изображения Landsat с помощью Python
Сопоставление положений пикселей с географическими координатами с использованием метаданных Landsat и Rasteriotowardsdatascience.com

Круг интересов

Во-первых, мы увидим, как закрасить AOI. Мы начинаем с создания маски, используя все 4 типа покрытия (строка 2), и удаляем ненужные пиксели с облачного изображения (строка 5). Затем мы обрезаем полученное изображение вокруг нашей области интереса (строки 8–9). Мы принимаем (x, y) за центр области интереса на облачном изображении.

# get mask for cloudy image
mask = cloud_mask*shadow_mask*dilated_cloud_mask*cirrus_mask

# remove cloudy pixels and fill with adjusted clear pixels
rm_mask = cloudy*mask[:, :, np.newaxis]

#crop area of interest
x,y,h = 4500,4500,300
crop_rm_mask = rm_mask[y-h:y+h,x-h:x+h]

Как уже упоминалось, нам нужно учитывать геолокации пикселей. Мы используем пакет rasterio, чтобы помочь с этим. Мы загружаем красную полосу как для облачного (строка 4), так и для ясного (строка 5) изображения.

import rasterio as rio

# Open red band with rasterio for geolocation
geo_cloudy = rio.open('./data/{}/{}_SR_B4.TIF'.format(IDs[0], IDs[0]))
geo_clear = rio.open('./data/{}/{}_SR_B4.TIF'.format(IDs[1], IDs[1]))

Мы находим UTM-координаты на облачном изображении центра нашей области интереса (строка 2). Затем мы находим позиции массива для этих координат UTM в чистом изображении (строка 5). Их печать дает значение (4430, 4510). Вы можете видеть, как это отличается от исходного центра (4500,4500).

# get UTM coordinates from cloudy image
utmx, utmy = geo_cloudy.xy(y,x)

# get pixels from clear image
y_,x_ = geo_clear.index(utmx, utmy)
print(x_,y_)

Мы вырезаем AOI из чистого изображения (строка 2), используя новый центр. Затем мы обрезаем облачную маску (из облачного изображения), используя исходный центр (строка 5). Координаты UTM crop_clear и crop_mask будут идеально совмещены. Итак, мы создаем нашу заливку, умножая чистое изображение на инверсию маски (строка 6). Последний шаг — добавить заливку к облачному изображению без облаков (строка 9).

# crop clear image 
crop_clear = clear[y_-h:y_+h,x_-h:x_+h]

#get fill from clear image
crop_mask = mask[y-h:y+h,x-h:x+h]
fill = crop_clear*(1-crop_mask[:, :, np.newaxis])

#inpaint area of interest 
inpaint = crop_rm_mask + fill

Результат вы можете увидеть на Рис. 8. Выглядит хорошо! По большей части пейзаж выглядит сплошным. Хотя есть несколько вопросов, которые мы обсудим в конце.

Все изображение

Закрашивать всю сцену сложнее. Для этого нам нужно выровнять ограничивающую рамку облачной и ясной сцен. Связанный прямоугольник — это внешняя черная граница сцены. Мы выравниваем UTM-координаты блоков, добавляя и удаляя некоторые пиксели из чистого изображения.

Начнем с расчета количества пикселей, которые нужно добавить/удалить сверху, снизу, слева и справа от четкого изображения. Мы получаем позиции массива для верхнего левого (UL) (строка 3) и нижнего правого (LR) углов (строка 4) чистого изображения. Они основаны на размерах четкого изображения (строка 2).

# get pixel coordinates of clear image corners
y,x = geo_clear.read(1).shape
clear_ul = (0,0) # upper left
clear_lr = (y,x) # lower right

Затем мы получаем границы из облачного изображения (строка 2). Они дают UTM-координаты ограничивающей рамки облачного изображения. Затем мы получаем позиции массива для этих координат из чистого изображения (строки 3–4). Обратите внимание, что эти позиции могут выходить за пределы четкого изображения.

# get pixel coordinates of cloudy image corners
cloudy_bounds = geo_cloudy.bounds
new_clear_ul = geo_clear.index(cloudy_bounds.left,cloudy_bounds.top)
new_clear_lr = geo_clear.index(cloudy_bounds.right,cloudy_bounds.bottom)

Разница между исходным и новым угловым положением подскажет нам, как нам нужно настроить четкое изображение, чтобы выровнять координаты UTM. Печать настроек (строка 8) дает нам -10 0 70 -70. Это означает, что мы должны настроить четкое изображение:

Удаление 10 пикселей сверху
Нет регулировки по низу
Добавление 70 пикселей влево
Удаление 70 пикселей справа

# calculate pixel adjustment 
top_adj = clear_ul[0] - new_clear_ul[0]
bottom_adj = new_clear_lr[0] - clear_lr[0]

left_adj = clear_ul[1] -  new_clear_ul[1]
right_adj = new_clear_lr[1] - clear_lr[1]

print(top_adj, bottom_adj, left_adj, right_adj)

После расчета корректировок мы можем использовать функцию adjust_rgb. Это добавит черные пиксели в случае положительных настроек (строки 6–17). Затем он обрежет пиксели для негативных настроек (строки 20–27).

def adjust_rgb(rgb,top_adj, bottom_adj, left_adj, right_adj):

    adj_rgb = rgb.copy()

    #Adding black pixels 
    if top_adj > 0:
        add_top = np.zeros((top_adj,rgb.shape[1],3))
        adj_rgb = np.vstack((add_top,adj_rgb))
    if bottom_adj > 0:
        add_bottom = np.zeros((bottom_adj,rgb.shape[1],3))
        adj_rgb = np.vstack((adj_rgb,add_bottom))
    if left_adj > 0:
        add_left = np.zeros((rgb.shape[0],left_adj,3))
        adj_rgb = np.hstack((add_left,adj_rgb))
    if right_adj > 0:
        add_right = np.zeros((rgb.shape[0],right_adj,3))
        adj_rgb = np.hstack((adj_rgb,add_right))
    
    #Removing pixels
    if top_adj < 0:
        adj_rgb = adj_rgb[-top_adj:,:,:]
    if bottom_adj < 0:
        adj_rgb = adj_rgb[:bottom_adj,:,:]
    if left_adj < 0:
        adj_rgb = adj_rgb[:,-left_adj:,:]
    if right_adj < 0:
        adj_rgb = adj_rgb[:,:right_adj,:]

    return adj_rgb

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

# Get RGB images
cloudy_RGB = get_RGB(IDs[0])
clear_RGB = get_RGB(IDs[1])

# Adjust clear RGB image
clear_RGB_adj = adjust_rgb(clear_RGB,top_adj, bottom_adj, left_adj, right_adj)

# get mask for cloudy image
mask = cloud_mask*shadow_mask*dilated_cloud_mask*cirrus_mask

# remove cloudy pixels and fill with adjusted clear pixels
rm_mask = cloudy_RGB*mask[:, :, np.newaxis]
fill_mask = clear_RGB_adj*(1-mask[:, :, np.newaxis])
inpaint = rm_mask+fill_mask

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

Мы можем увидеть это, если увеличим масштаб другого AOI. Прозрачные белые пиксели — это те, которые были заменены пикселями в чистом изображении. Обратите внимание, что некоторые теневые пиксели не были заменены. Это может привести к разрывам в пейзаже, где резко меняется яркость.

Проблема в том, что облачные маски не идеальны. Мы не можем полностью заменить все облака, если не можем идентифицировать их со 100% точностью. Разрывы также могут быть вызваны изменениями ландшафта с течением времени. Реки движутся, и приливы будут расти. Чтобы создать наилучшую визуализацию, выберите заменяющие изображения, которые были сделаны близко ко времени вашего облачного изображения.

Надеюсь, вам понравилась эта статья! Вы можете найти меня на Mastodon | Твиттер | Ютуб | Информационный бюллетень — подпишитесь на БЕСПЛАТНЫЙ доступ к Курсу Python SHAP

Присоединяйтесь к Medium по моей реферальной ссылке — Конор О'Салливан
Как участник Medium, часть вашего членского взноса идет на авторов, которых вы читаете, и вы получаете полный доступ ко всем история…conorosullyds.medium.com