-данные Sentinel NDVI в качестве примера
TL;DR
- Добавьте свой шейп-файл в активы
- Используйте API JavaScript Google Earth Engine.
- Создайте нормализованный разностный индекс растительности (NDVI) на основе данных Sentinel.
- Извлечение и визуализация временных рядов в консоли редактора кода Earth Engine
- Извлечение и экспорт больших данных временного ряда в CSV
При работе с геопространственными данными нас часто интересует извлечение временных рядов для набора местоположений. Это, например, могут быть несколько точек или многоугольников, представляющих разные местоположения. В этом коротком блоге мы будем использовать JavaScript в редакторе кода Google Earth Engine и извлекать временные ряды по набору полигонов. Мы будем использовать данные Sentinel-2 для извлечения индекса нормализованной разности растительности (NDVI) для набора орошаемых полей. Давайте начнем с переноса наших полигональных данных в движок Google Earth.
Добавление шейп-файла многоугольника к ресурсам
Чтобы загрузить шейп-файл в качестве ресурса в Google Earth Engine, перейдите к Assets на левой боковой панели, нажмите New и выберите Shapefiles в раскрывающемся списке. Появится окно для загрузки шейп-файла, как показано на снимке экрана ниже.
Select файл с вашего компьютера для загрузки и нажмите загрузить. В этом случае файл называется mgck_shp.zip. Вы сможете увидеть ход загрузки на боковой панели Tasks справа (см. снимок экрана ниже), а когда загрузка будет завершена, шейп-файл появится под CLOUD ASSETS на левой боковой панели (в данном случае выделен желтым цветом).
Импорт и визуализация FeatureCollection
Только что добавленный в Assets шейп-файл можно загрузить и отобразить на карте с помощью следующего фрагмента кода:
//Import GEE Feature Collection (field geometry)
var geometry =ee.FeatureCollection("projects/ee-chinmaydeval91/assets/mgck_shp");
//Display polygons on a map
Map.addLayer(geometry); Map.centerObject(geometry,14);
На карте показаны восемь полей орошаемого земледелия в Мэджик-Вэлли, штат Айдахо.
Создание и извлечение Sentinel NDVI
Насколько мне известно, в отличие от составных продуктов MODIS NDVI, для Sentinel на GEE нет такого предварительно вычисленного продукта NDVI. Поэтому сначала нам нужно создать NDVI.
Импорт коллекции изображений
Начнем с импорта коллекции изображений отражательной способности поверхности Sentinel-2, введя следующий код:
// Import image collection of Sentinel-2 imagery, and
// filter by start and end date, and
// filter by boundary extent
var S2 = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2019-03-28', '2021-12-31')
.filterBounds(geometry);
Это импортирует коллекции Sentinel-2 и фильтрует их по заданному диапазону дат, а также по границам коллекции объектов. Далее нам нужно замаскировать облачные и снежные изображения.
Маскируйте облака и снег
Я наткнулся на эту функцию maskS2sr в Руководстве по составлению карт UN Spider’s Burn Severity, которая написана для маскировки облаков и снега. Я буду использовать его здесь:
// Function to mask clouds from the pixel quality band of Sentinel-2 SR data.
function maskS2sr(image) {
// Bits 10 and 11 are clouds and cirrus, respectively.
var cloudBitMask = ee.Number(2).pow(10).int();
var cirrusBitMask = ee.Number(2).pow(11).int();
// Get the pixel QA band.
var qa = image.select('QA60');
// All flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
.and(qa.bitwiseAnd(cirrusBitMask).eq(0));
// Return the masked image, scaled to TOA reflectance, without the QA bands.
return image.updateMask(mask)
.copyProperties(image, ["system:time_start"]);
}
// mask clouds and snow across the image collection
var S2 = S2.map(maskS2sr)
Теперь у нас есть коллекция изображений Sentinel-2, свободная от облаков и снега и отфильтрованная для заданных пространственных и временных границ. Мы даже можем прикрепить коллекцию изображений к этим орошаемым полям.
Привязка к полигонам орошаемых полей
Чтобы обрезать изображение по полигону, мы можем использовать следующую функцию clipImgCollect, а затем map по всей коллекции изображений:
////// function to clip the entire ImageCollection to the field boundary //////
function clipImgCollect(img) {
return img.clip(geometry);
}
// clip S2 ImageCollection to poly bounds
var clippedS2 = S2.map(clipImgCollect);
Рассчитать NDVI
Теперь последний шаг перед визуализацией — расчет NDVI. Это можно сделать с помощью приведенной ниже функции NDVI_from_S2, которая выполняет математические вычисления диапазона. Затем мы можем map использовать эту функцию для всей коллекции изображений:
// Function to calculate and create NDVI
var NDVI_from_S2 = function(img) {
return img.addBands(img.normalizedDifference(['B8', 'B4']));
};
// map function across the S2 image collection and
// create NDVI
var clippedS2 = clippedS2.map(NDVI_from_S2);
Этот шаг добавит полосу nd, содержащую NDVI, в коллекцию clippedS2.
Визуализируйте временные ряды
Чтобы создать график временных рядов, используйте функцию seriesByRegion и передайте статистику для расчета методу ee.Reducer, который помогает агрегировать данные. В этом случае мы будем агрегировать средние значения данных по геометрии с помощью функции ee.Reducer.mean().
// Create a time series chart.
var plotNDVI_TS = ui.Chart.image.seriesByRegion(
clippedS2,
geometry,ee.Reducer.mean(),
'nd',
500,
'system:time_start', 'system:index')
.setChartType('LineChart')
.setOptions({ // Plot customization options
interpolateNulls: true,
lineWidth: 1,
pointSize: 2,
title: 'Field mean NDVI',
hAxis: {title: 'Date'},
vAxis: {title: 'NDVI'}
});
print(plotNDVI_TS)
Это создаст график среднего полевого временного ряда NDVI для всех орошаемых полей, как показано ниже. Эту диаграмму можно загрузить как есть или даже загрузить данные, нанесенные на диаграмму, в виде таблицы (CSV).
Однако создание plotNDVI_TS, как показано выше, безусловно, завершится ошибкой, если ваш запрос запроса превышает 5000 элементов. Например, в этом примере я просто увеличил размер своей коллекции изображений, выбрав более старую дату начала (2017-03-28) в функции filterDate. Когда я выполнил функцию plotNDVI_TS, я получил в вашей консоли ошибку, которая выглядит так:
Запрашиваемые данные слишком велики для отображения в консоли, и их необходимо экспортировать из механизма Google Earth.
Экспорт временных рядов в табличном формате
Чтобы иметь возможность экспортировать большие временные ряды, сначала выберите полосу NDVI из коллекции изображений clippedS2. Затем сопоставьте функцию reduceRegions с этой коллекцией изображений NDVI, чтобы вычислить среднее значение NDVI для каждого объекта (полигона орошаемого поля) для всех изображений в коллекции. Мы также сохраним дату изображения в свойстве date и отформатируем ее как ГГГГ-ММ-ДД. Этот объект dat даст нам идентификатор полигона, дату изображения и среднее значение NDVI для этого полигона.
// Select ndvi band
var clippedS2_ndvi = clippedS2.select(['nd'])
// Collect image date and ndvi mean value for all features
var dat = clippedS2_ndvi.map(function(image) {
return image.reduceRegions({
collection:geometry ,
reducer: ee.Reducer.mean().setOutputs(['ndvi']),
scale: 10
}).filter(ee.Filter.neq('ndvi', null))
.map(function(f) {
return f.set('date', image.date().format('YYYY-MM-dd'));
});
}).flatten();
Теперь, чтобы подготовить таблицу к экспорту, используйте функцию format, показанную ниже, и примените ее к объекту dat, чтобы отформатировать его таким образом, чтобы каждая строка в таблице содержала временной ряд среднего значения NDVI для каждого полигона орошаемого поля. Заголовки столбцов для каждого из этих значений представляют собой даты изображений, для которых было рассчитано среднее значение NDVI. Отформатированная таблица хранится в объекте sentinelNDVI_Results. Эта функция format была адаптирована из источника, связанного здесь.
var format = function(table, rowId, colId) {
var rows = table.distinct(rowId);
var joined = ee.Join.saveAll('matches').apply({
primary: rows,
secondary: table,
condition: ee.Filter.equals({
leftField: rowId,
rightField: rowId
})
});
return joined.map(function(row) {
var values = ee.List(row.get('matches'))
.map(function(feature) {
feature = ee.Feature(feature);
var ndvi_val = ee.List([feature.get('ndvi'), -9999]).reduce(ee.Reducer.firstNonNull());
return [feature.get(colId), ee.Number(ndvi_val).format("%.3f")];
});
return row.select([rowId]).set(ee.Dictionary(values.flatten()));
});
};
// Apply format function to the dat object
var sentinelNDVI_Results = format(dat, 'ID', 'date');
Экспорт этой отформатированной таблицы на Google Диск довольно прост. Просто передайте объект sentinelNDVI_Results функции Export.table.toDrive. Другие аргументы, которые вам необходимо указать, включают folder на диске, на который вы хотите экспортировать файл, fileNamePrefix — префикс, который будет использоваться в имени файла, и fileFormat:
// Export data frame to Google Drive as a CSV
Export.table.toDrive({
collection: sentinelNDVI_Results,
description: 'Sentinel_NDVI_TS',
folder: 'earthengine',
fileNamePrefix: 'sentinel_ndvi_ts_demo_',
fileFormat: 'CSV'
})
Вот краткий обзор экспортированного временного ряда, который я быстро построил в Excel для демонстрации:
Вы можете найти полный код, используемый в этом блоге, по ссылке эта ссылка GEE.
Первоначально опубликовано на https://chinmaydeval.com.
