Для получения более подробной информации вы можете посмотреть видео на YouTube о том, как вычислить статистику спутниковых изображений с помощью pandas в google colab.
URL-адрес видео:https://youtu.be/xPTAl_qf864
Обзор
В этом руководстве мы рассмотрим, как вычислить статистику спутниковых изображений с помощью pandas в google colab и получить зональную статистику из изображений Sentinel 2, которые мы объединим с нашим фреймом данных pandas, в пошаговом руководстве о том, как рассчитать водный индекс NDWI. для затопленных территорий.
Спутниковые данные являются плотными при использовании в искусственном интеллекте машинного обучения, и данные хранятся в виде значений с использованием ячеек. Но нам нужна только сводка спутникового изображения, известного как растровое изображение, которое будет преобразовано в табличный формат в CSV или Pandas. кадр данных. Спутниковые изображения представляют собой попиксельные данные, называемые Растрами.
Растровые изображения в основном состоят из спутниковых изображений и растров, состоящих из матрицы ячеек и строк. Каждая ячейка или строка содержит информацию о местоположении, такую как высота, температура и растительность. Этот процесс получения табличных результатов, называемый Сводная статистика из растровых изображений, называется Зональной статистикой.
Исследование данных
Во-первых, установите и импортируйте библиотеки.
!apt install gdal-bin python-gdal python3-gdal !pip install rasterio !apt install python3-rtree !pip install git+git://github.com/geopandas/geopandas.git !pip install descartes !pip install mapclassify !pip install sentinelsat import pandas as pd import numpy as np import geopandas as gpd import rasterio as rio from rasterio.plot import show from rasterio.mask import mask import matplotlib.pyplot as plt import seaborn as sns
Gdal — популярная и мощная библиотека обработки геопространственных данных, доступная в виде набора утилит командной строки (готовых к использованию двоичных файлов).
Mapclassify – это библиотека python с открытым исходным кодом для классификации карт Choropleth. Это часть рефакторинга PySAL.
Sentinelsat — библиотека Python, упрощающая поиск, извлечение и загрузку спутниковых изображений Sentinel. Итак, давайте начнем установку sentinelsat через pip.
Итак, далее мы будем загружать файлы из Интернета и распаковывать их.
!wget https://www.dropbox.com/s/gf51dybqbujyjb2/Data.zip !unzip Data.zip
Теперь с помощью библиотеки Rasterio мы можем читать различные полосы изображений Sentinel 2. В этом случае мы читаем полосы 8 как NIR, 4 как красные, 3 как зеленые и 2 как синие.
b8 = rio.open("/content/Data/20191101/B08-20191101.tif")
b4 = rio.open("/content/Data/20191101/B04-20191101.tif")
b3 = rio.open("/content/Data/20191101/B03-20191101.tif")
b2 = rio.open("/content/Data/20191101/B02-20191101.tif")
Далее Давайте посмотрим на ширину и высоту изображений. Сейчас мы будем использовать только b4, но мы можем проверить, все ли ленты имеют одинаковый вес и длину или нет.
b4.count, b4.width, b4.height
Мы можем построить данные, чтобы увидеть и изучить спутниковые изображения, которые у нас есть. Здесь мы будем строить график только для Band 3.
fig, ax = plt.subplots(1, figsize=(12, 10)) show(b3, ax=ax) plt.show()
Изображение области, сделанное Sentinel 2 только для диапазона 3, показано ниже.
Давайте также прочитаем таблицу зданий, которую мы будем использовать для хранения статистических сводок, полученных из спутникового снимка.
Обратите внимание, мы можем использовать другие полигоны, такие как районы, прямоугольные сетки вместо полигонов зданий. Мы будем использовать Geopandas для чтения данных.
buildings = gpd.read_file("/content/Data/shapefiles/osm_buildings.shp")
buildings = buildings[["osm_id", "building", "geometry"]]
buildings.head()
А вот и первые пять строк таблицы, показанной ниже. У нас есть столбец геометрии каждого здания, osm_id и столбцы здания.
Мы также можем нанести здания поверх изображений Sentinel 2.
fig, ax = plt.subplots(figsize=(12, 10)) show(b4, ax=ax) buildings.plot(ax=ax, color="white", alpha=.50) plt.show();
Здания отмечены белым цветом и перекрываются слоями на изображении, как показано на рисунке. На изображении показаны размеры города, т. е. жилые районы с извилистой рекой Шабелле.
Расчет NDWI (Рассчитать нормализованный разностный индекс воды)
Чтобы рассчитать значения NDWI по изображениям Sentinel 2, мы используем формулу:
(Band3 — Band8)/(Band3 + Band8)
Итак, давайте посчитаем по этой формуле в Rasterio.
green = b3.read() nir = b8.read() ndwi = (nir.astype(float)-green.astype(float))/(nir+green)
Массивы NDWI можно построить с помощью Rasterio, как показано ниже.
fig, ax = plt.subplots(1, figsize=(12, 10)) show(ndwi, ax=ax, cmap="coolwarm_r") plt.show()
На графике NDWI четко показаны затопленные районы у реки Шабель синим цветом, а степень затопления приходится на жилые районы.
Мы можем сохранить массивы NDWI как растровое изображение, чтобы использовать его позже.
meta = b4.meta
meta.update(driver='GTiff')
meta.update(dtype=rio.float32)
with rio.open('NDWI.tif', 'w', **meta) as dst:
dst.write(ndwi.astype(rio.float32))
ndwi_raster = rio.open('NDWI.tif')
Теперь, когда мы рассчитали значения NDWI, пришло время получить статистику из растрового изображения NDWI и объединить ее с нашей таблицей зданий. Мы используем функции маски Rasterio, чтобы получить значения ячеек из растрового изображения NDWI.
Ниже приведена небольшая функция, которая маскирует значения ячеек для нашей таблицы фрейма данных.
def derive_stats(geom, data=ndwi_raster, **mask_kw): masked, mask_transform = mask(dataset=data, shapes=(geom,), crop=True, all_touched=True, filled=True) return masked
Теперь мы можем получить нужные нам значения. Допустим, нам нужно получить средние значения NDWI для каждого здания. Мы создаем столбец для этого «mean_ndwi» и передаем нашу функцию для применения геометрии здания, а также применяем среднее значение с помощью Numpy.
buildings['mean_ndwi'] = buildings.geometry.apply(derive_stats).apply(np.mean)
Или получите максимальные значения NDWI для каждого здания.
buildings['max_ndwi'] = buildings.geometry.apply(derive_stats).apply(np.max) buildings.head()
В нашей таблице есть два новых столбца, т. е. mean_ndwi и max_ndwi, где мы можем хранить средние и максимальные значения NDWI для каждого здания.
Давайте создадим фигуру и набор подграфиков для объекта оси, показав среднее значение NDWI.
fig, ax = plt.subplots(figsize=(12,10)) buildings.plot(column="mean_ndwi", ax=ax, cmap="Reds_r", scheme='quantiles', legend=True)
Затем давайте создадим график морского побережья для максимального NDWI.
sns.distplot(buildings.max_ndwi)
Затем начертите фигуру max_ndwi через растр ndwi жилых районов темно-зеленым цветом.
fig, ax = plt.subplots(figsize=(12, 10)) show(ndwi_raster, ax=ax) buildings[buildings["max_ndwi"].between(0.4, 0.6)].plot(color= "darkgreen", ax=ax) plt.show();
Наконец, мы построим фигуру max_ndwi через растр ndwi жилых районов красным цветом.
fig, ax = plt.subplots(figsize=(12, 10)) show(ndwi_raster, ax=ax) buildings[buildings["max_ndwi"].between(0, 0.1)].plot(color="red", ax=ax) plt.show();
