Упростите загрузку сцен Landsat с помощью пакета Python landsatxplore
Спутники Landsat являются одними из наиболее часто используемых источников данных наблюдения Земли. Они предоставляют высококачественные изображения поверхности планеты уже более четырех десятилетий. Однако ручная загрузка этих изображений может быть утомительной! К счастью, с пакетом landsatxplore вы можете легко загружать и обрабатывать сцены Landsat с помощью нескольких строк кода.
Мы изучим пакет landsatxplore и пройдем шаги по загрузке спутниковых изображений Landsat с помощью Python. Это включает в себя:
- Настройка доступа к API с учетной записью USGS
- Поиск и фильтрация сцен Landsat
- Загрузка и работа со сценами в Python
Попрощайтесь с ручными загрузками и поприветствуйте автоматизированный и эффективный рабочий процесс!
Настройка landsatexplore
Шаг 1: Зарегистрируйтесь в Геологической службе США
Для начала вам нужно будет настроить учетную запись USGS. Это та же учетная запись, которую вы используете для загрузки сцен с EarthExplorer. Сохраните свое имя пользователя и пароль, так как они понадобятся нам позже.
После регистрации вы можете использовать USGS M2M API. Однако для его настройки потребуется некоторая работа. Вместо этого мы будем использовать пакет lansatxplore. Он абстрагирует большую часть технических деталей для вас.
Шаг 2: Установите lansatxplore
Следуйте инструкциям на GitRepo.
Шаг 3. Проверьте соединение API
Используйте приведенный ниже код, чтобы убедиться, что все настроено правильно. Вы должны заменить имя пользователя и пароль данными, которые вы использовали для регистрации своей учетной записи Геологической службы США.
from landsatxplore.api import API # Your USGS credentials username = "XXXXXXXXXXXX" password = "XXXXXXXXXXXX" # Initialize a new API instance api = API(username, password) # Perform a request response = api.request(endpoint="dataset-catalogs") print(response)
Вывод ответа должен выглядеть следующим образом:
{‘EE’: ‘EarthExplorer’, ‘GV’: ‘GloVis’, ‘HDDS’: ‘HDDS Explorer’}
Это наборы данных, доступные через API. В этом руководстве мы рассматриваем только набор данных EarthExplorer.
Поиск по сценам с помощью EarthExplorer
Прежде чем мы перейдем к загрузке сцен с помощью API, мы выполним ручной поиск с помощью EarthExplorer. Это сделано для того, чтобы мы могли сравнить результаты с тем, что мы видим при использовании Python. Если вы не знакомы с порталом EarthExplorer, вам может помочь это руководство.
Мы ищем сцены по следующим критериям:
- Они должны содержать точку на заданной широте и долготе. Эта точка приходится на остров Бык в Дублине.
- Снято с 01.01.2020 по 31.12.2022
- Максимальная облачность 50%
- Часть коллекции Landsat 8 или 9 уровня 2.
Результаты поиска вы можете увидеть ниже. Мы отмечаем некоторые вещи для сравнения с нашим поиском Python:
- Есть 54 сцены, соответствующие критериям поиска.
- На острове Бык есть 2 плитки с точкой. Они имеют значения пути и строки (206, 023) и (205 023).
- Идентификатор первой сцены: LC08_L2SP_206023_20221118_20221128_02_T1. Если вам интересно, что означает этот идентификатор, ознакомьтесь с Соглашением об именах Landsat.
Использование пакета lansatxplore Python
Поиск сцен
Выполним аналогичный поиск с помощью landsatxplore. Мы делаем это, используя функцию поиска ниже:
- набор данных — определяет, для какого спутника нам нужны сцены. Мы использовали идентификатор набора данных для Landsat 8 и 9. См. GitRepo для идентификатора для Landsat 5 и 7.
- широта и долгота дают одну и ту же точку на острове Бык. За исключением того, что мы перевели координаты в десятичные дроби.
- start_date, end_date и max_cloud_cover также остались прежними.
# Search for Landsat TM scenes
scenes = api.search(
dataset='landsat_ot_c2_l2',
latitude=53.36305556,
longitude=-6.15583333,
start_date='2020-01-01',
end_date='2022-12-31',
max_cloud_cover=50
)
# log out
api.logout()
Результат поиска вернет информацию в словаре JSON. Мы конвертируем это в Pandas DataFrame (строка 4), где каждая строка будет представлять уникальную сцену. Возвращается много метаданных! Итак, фильтруем то, что необходимо для данного приложения (строка 5). Наконец, мы упорядочены по acquisition_date — дате съемки сцены спутником Landsat.
import pandas as pd
# Create a DataFrame from the scenes
df_scenes = pd.DataFrame(scenes)
df_scenes = df_scenes[['display_id','wrs_path', 'wrs_row','satellite','cloud_cover','acquisition_date']]
df_scenes.sort_values('acquisition_date', ascending=False, inplace=True)
Вы можете увидеть снимок этого набора данных ниже. Сравнивая это с нашим поиском с помощью EarthExplorer, мы можем быть уверены, что результаты такие же. Этот набор данных содержит 54 строки и две уникальные пары wrs_path и wrs_row — (206, 23) и (205,23).Первый display_id такой же, как мы видели раньше.
Если бы мы захотели, мы могли бы отфильтровать этот набор данных дальше. Мы могли бы использовать столбец satellite, чтобы выбрать только изображения со спутников Landsat 8 или 9. Кроме того, в столбце cloud_cover указан процент покрытия изображения облаками. Когда вы будете довольны окончательным списком сцен, вы можете перейти к их загрузке.
Загрузка данных
Ниже у нас есть код, используемый для загрузки сцены Landsat. Мы используем функцию EarthExplorer (строка 5). Это инициализируется так же, как и раньше — с использованием ваших учетных данных USGS. Чтобы загрузить сцену, нам нужно использовать ее display_id и определить выходной каталог (строка 12). Мы используем display_id первой сцены, упомянутой выше (строка 8).
from landsatxplore.earthexplorer import EarthExplorer
import os
# Initialize the API
ee = EarthExplorer(username, password)
# Select the first scene
ID = 'LC08_L2SP_206023_20221118_20221128_02_T1'
# Download the scene
try:
ee.download(ID, output_dir='./data')
print('{} succesful'.format(ID))
# Additional error handling
except:
if os.path.isfile('./data/{}.tar'.format(ID)):
print('{} error but file exists'.format(ID))
else:
print('{} error'.format(ID))
ee.logout()
Возможно, вы заметили дополнительную обработку ошибок выше. Это из-за проблемы с пакетом. В некоторых случаях сцены загружаются правильно, но ошибка все равно возникает. Дополнительная обработка ошибок дважды проверяет существование файла сцены.
Работа с данными
Сцена будет загружена в виде файла tar. Имя файла будет display_id, за которым следует .tar:
LC08_L2SP_206023_20221118_20221128_02_T1.tar
Мы можем работать с этими данными прямо в Python. Для начала нам нужно разархивировать tar-файл (строки 4–6). Имя новой папки устанавливается равным display_id сцены (строка 5).
import tarfile
# Extract files from tar archive
tar = tarfile.open('./data/{}.tar'.format(ID))
tar.extractall('./data/{}'.format(ID))
tar.close()
Вы можете увидеть разархивированную папку и все доступные файлы ниже. Сюда входит вся информация, доступная для Лансата научные продукты уровня 2. Применения для этих данных бесконечны! Например, мы визуализируем эту сцену, используя полосы видимого света. Они доступны в файлах, выделенных ниже.
Загружаем синюю, зеленую и красную полосы (строки 6–8). Складываем эти бэнды (строка 11), масштабируем (строка 12) и обрезаем для усиления контраста (строка 15). Наконец, мы отображаем это изображение с помощью matplotlib (строки 18–20). Вы можете увидеть это изображение ниже.
import tifffile as tiff
import numpy as np
import matplotlib.pyplot as plt
# 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.2)/0.2
# Display RGB image
fig, ax = plt.subplots(figsize=(10, 10))
plt.imshow(RGB)
ax.set_axis_off()
Прочтите эту статью, если хотите узнать больше о работе с RGB-каналами спутниковых изображений:
Визуализация каналов RGB — это только начало. После загрузки данных мы можем выполнять любые задачи ДЗЗ — от расчета индексов до обучения моделей. Самое приятное то, что мы можем сделать все это, не покидая блокнота Python.
Надеюсь, вам понравилась эта статья! Вы можете поддержать меня, став одним из моих приглашенных участников :)
| Твиттер | Ютуб | Информационный бюллетень — подпишитесь на БЕСПЛАТНЫЙ доступ к Курсу Python SHAP
