Более продвинутое использование куба данных

Построение одиночного диапазона данных

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

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

import datacube
import matplotlib.pyplot as plt

dc = datacube.Datacube()

query = {'x': [75.0,75.5],
         'y': [40.0,40.5],
         'time': ['2019-07-01','2019-07-10'],
         'measurements': ['red','green','blue','mask'],
         'output_crs': 'EPSG:32643',
         'resolution': [20,20]
        }

ds = dc.load(product='s2_10m', **query)
print(ds)

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

td = ds.isel(time=2)
print(td)
td.red.plot()
plt.show()

Надеюсь, что изображение было показано выше с использованием matplotlib для его отображения (модуль matplotlib также был рассмотрен в третьем блокноте этой серии). Теперь я опишу команду для отображения красной полосы, ds.red.plot(). ds сам по себе не может быть нанесен на график, так как это объект набора данных и потенциально может содержать данные за несколько дней или несколько каналов. Однако, чтобы выбрать только одну полосу в качестве DataArray, представляющего собой объект, который можно построить, вы должны ввести команду ds."band_desired_to_be_selected". Поэтому, если бы я просто хотел выбрать, например, зеленую полосу, я бы ввел ds.green.

Любой объект DataArray, который представляет собой всего лишь один временной шаг, теоретически может быть построен как изображение с помощью метода .plot(). Это также будет работать, если у вас уже есть объект DataArray, например, сохраненный как переменная. Также можно строить изображения RGB, используя видимые каналы или создавая составное изображение в ложных цветах. Это немного более сложная задача, и мы ее сейчас рассмотрим.

Создание RGB изображения

Создание RGB изображения требует немного больше работы, чем просто вызов .plot () для DataArray, поскольку теперь необходимо построить три полосы. Объект ds, который мы построили выше, содержит четыре канала, поскольку маска облака присутствует в наборе данных, поэтому мы изменим запрос, чтобы получить только три канала из куба данных, выполнив приведенный ниже код:

query = {'x': [75.0,75.5],
         'y': [40.0,40.5],
         'time': ['2019-07-01','2019-07-10'],
         'measurements': ['red','green','blue'],
         'output_crs': 'EPSG:32643',
         'resolution': [20,20]
        }

ds = dc.load(product='s2_10m', **query)
print(ds)

Первое, что нужно сделать для построения изображения RGB, - это превратить Dataset в DataArray с одним временным шагом, но тот, который по-прежнему содержит три требуемые полосы. Для этого мы сначала выбираем один временной шаг с помощью метода .isel(). Это требуется даже в приведенном ниже примере, где есть только один временной шаг в наборе данных, поскольку ..isel() полностью удаляет его как измерение массива. Затем массив можно преобразовать в DataArray с помощью команды .to_array() и указания измерения, по которому теперь будут разделяться исходные полосы в наборе данных - в этом случае мы называем это измерение цветовым измерением.

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

Эти описанные шаги показаны ниже:

ds = ds.isel(time=2)
rgb = ds.to_array(dim='color')

fake_saturation = 6000

Следующие несколько шагов довольно сложны. Из-за того, как устроен DataArray, на данный момент размеры изображения упорядочены ("color", "x", "y"). Чтобы построить изображение RGB с помощью Python, необходимо указать измерение «цвет» после размеров «y» и «x» в вызове функции plot. Поэтому нам нужно переместить его, и это можно сделать с помощью метода .transpose(). Команда, показанная ниже, изменит порядок размеров с ("color", "x", "y") на ("y","x","color"):

rgb = rgb.transpose(*(rgb.dims[1:]+rgb.dims[:1]))

Следующие два шага описывают масштабирование данных, чтобы убедиться, что ничто не превышает заданную нами переменную fake_saturation, а также для обеспечения масштабирования всех данных в полосах изображения от 0 до 1, поскольку это требуется для построения изображения в формате RGB. Чтобы гарантировать, что ни один пиксель не превышает уровень fake_saturation, мы маскируем данные, чтобы любые пиксели выше этого значения были равны значению. Это делается с помощью метода .where(*insert condition*).

Затем все пиксели в DataArray делятся на значение fake_saturation. Поскольку максимальное значение массива теперь равно fake_saturation, а минимальное значение равно 0, это означает, что все значения в DataArray теперь находятся между 0 и 1. Это делается с помощью оператора /=, который делит значение слева. рядом со значением в правой части уравнения и устанавливает результат в качестве нового значения в левой части.

Эти два шага показаны ниже:

rgb = rgb.where((rgb <= fake_saturation).all(dim='color'))
rgb /= fake_saturation
print(rgb)

Теперь у нас есть DataArray, который готов к построению в виде изображения RGB. Это можно сделать аналогично тому, как это было сделано выше, за исключением того, что в этом случае мы используем метод .plot.imshow() вместо просто .plot(). Это связано с тем, что imshow поддерживает сбор и построение многополосных изображений.

rgb.plot.imshow()
plt.show()

В качестве альтернативы вы можете использовать ключевое слово robust=True в imshow для отображения данных. Это удалит аномальные пиксели, данные будут обрезаны на 2 и 98 процентилях.

ds[['red','green','blue']].to_array().plot.imshow(robust=True)
plt.show()

Создание индексов по спутниковым снимкам

Также возможно создание индексов по спутниковым снимкам. Куб данных Кыргызстана уже содержит некоторые средние медианные многодневные индексы, к которым можно получить доступ, изменив тип продукта на «индексы» при запросе куба данных.

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

query = {'x': [75.0,75.5],
         'y': [40.0,40.5],
         'time': ['2019-07-01','2019-07-10'],
         'measurements': ['red','nir','mask'],
         'output_crs': 'EPSG:32643',
         'resolution': [20,20]
        }

ds = dc.load(product='s2_10m', **query)
ds = ds.isel(time=2)
print(ds)

Теперь создать полосу NDVI из этих данных - несложная задача. NDVI рассчитывается с использованием следующего уравнения: $NDVI = \frac{NIR - RED}{NIR + RED}$

Используя полученные выше знания о том, что вы можете получить бэнды из набора данных с помощью ds."band_desired_to_be_selected", можно выполнить этот расчет для всех пикселей в бэнде следующим образом:

ndvi = ((ds.nir - ds.red) / (ds.nir + ds.red)).astype(float)
print(ndvi)

Метод .astype(float) в конце команды гарантирует, что выходной массив NDVI имеет тип float, так как странные результаты могут иногда возникать с Xarray при выполнении математических вычислений полного массива с целочисленными массивами.

Теперь массив NDVI можно построить как массив с одной полосой, вызвав для него метод .plot(), как это делалось ранее. Запустите приведенный ниже код, чтобы увидеть график NDVI.

ndvi.plot(cmap='RdYlGn')
plt.show()

В полосе маски каждый пиксель имеет значение, соответствующее разному типу значения маски. Например, чистая земля обозначается цифрой 5, снег / лед - 7, облака - 15 и т. Д. Полный список категорий можно найти в папке shared_data.

Что касается NDVI, нас интересуют только четкие пиксели земли, поэтому мы можем использовать один из методов, рассмотренных ранее .where(), чтобы замаскировать все, кроме земли. Как только это будет сделано, NDVI можно пересчитать и построить график, как показано ниже:

masked = ds.where(ds.mask == 5)
print(masked)

ndvi_masked = ((masked.nir - masked.red) / (masked.nir + masked.red)).astype(float)
ndvi_masked.plot(vmin=0,vmax=0.5,cmap='RdYlGn')
plt.show()