GEE-35. Временные ряды

Вопроса построения временных рядов мы уже касались. Речь шла о временных рядах для нескольких участков, но возможности GEE позволяют строить временные ряды одновременно для тысяч участков. При этом возникает несколько проблем, решением которых мы сейчас и займемся.

Район интереса geometry представлен тремя сельскохозяйственными полями 1️⃣. Для этих полей построим временные ряды нормализованного разностного вегетационного индекса (NDVI) за 2019 год — для каждого поля в каждую дату съемки рассчитывается среднее по участку значение NDVI. На практике такая задача возникает при мониторинге состояния сельскохозяйственных культур.

Сформируем коллекцию снимков Sentinel-2:

var col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
    .filterDate(startDate, endDate)
    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', cloudMax))
    .map(maskCloudAndShadows)
    .map(addNDVI)
    .filterBounds(geometry);

Здесь startDate и endDate — начальная и конечная даты периода наблюдений, cloudMax — максимальная доля облачного покрытия снимков. Функция maskCloudAndShadows маскирует облачность и тени облаков, а addNDVI добавляет к снимкам слой с рассчитанными значениями NDVI.

Запишем код, отвечающий за построение временных рядов:

var timeSeries = col.select('ndvi').map(function(image) {
  return image.reduceRegions({
    collection: geometry, 
    reducer: ee.Reducer.mean().setOutputs(['ndvi']).combine({
      reducer2: ee.Reducer.count(),
      sharedInputs: true
    }), 
    scale: 10,
  }).map(function(feature) {
    var ndvi = ee.List([feature.get('ndvi'), 'NA'])
      .reduce(ee.Reducer.firstNonNull());
    var imageID = image.id();
    var area = feature.area({'maxError': 1});
    var date = ee.String(imageID).slice(0,8);
    var share = ee.Number(feature.get('count')).multiply(10).multiply(10).divide(area);
    return feature.set({'id': feature.id(),'ndvi': ndvi, 'imageID': imageID, 'date':date, 'share': share});
    });
}).flatten();

Разбором кода займемся в следующем посте.

Код примера

#GEE

GEE-35. Временные ряды
Продолжение

Рассмотрим упрощенный код для построения временных рядов, который будем постепенно наращивать:

var timeSeries = col.select('ndvi').map(function(image) {
  return image.reduceRegions({
    collection: geometry, 
    reducer: ee.Reducer.mean().setOutputs(['ndvi']), 
    scale: 10,
  }).map(function(feature) {
    var ndvi = ee.List([feature.get('ndvi'), 'NA'])
      .reduce(ee.Reducer.firstNonNull());
    return feature.set({'ndvi': ndvi});
    });
}).flatten();

Из каждого снимка коллекции выбираем слои с NDVI и вычисляем средние значения индекса по каждому участку. В первой анонимной функции этим занимается reduceRegions. В результате получим ee.FeatureCollection, элементы которой содержат свойство 'ndvi'. Если бы мы не применили .setOutputs(['ndvi']), это свойство, по умолчанию, называлось бы mean.

К элементам полученной FeatureCollection применяется еще одна анонимная функция. Она нужна, чтобы свойство 'ndvi' всегда имело какое-то значение (кроме null). Дело в том, что reduceRegion не возвращает результат, если снимок не пересекается с полигоном района интереса или если фрагмент снимка удален после маскирования облачности. Для обработки таких ситуаций мы создаем список (ee.List), составленный из результата выполнения reduceRegion и пустого значения (‘NA’ — Not Available). Редьюсер ee.Reducer.firstNonNull() возвращает первое непустое значение. Если 'ndvi' окажется пустым (null), то этот редьюсер вернет ‘ndvi = NA’.

#GEE

GEE-35. Временные ряды
Продолжение

В отличие от данных MODIS и Landsat, гранулы Sentinel-2 пересекаются 1️⃣. Если район интереса имеет большую площадь, может возникнуть ситуация, когда один участок попадает на несколько соседних гранул Sentinel-2. То есть у одного участка в одну и ту же дату может оказаться несколько средних значений NDVI.

Решить эту проблему несложно. В качестве окончательного результата можно взять среднее или максимальное значение средних NDVI. Но для этого в готовой FeatureCollection нужно возвращать не только NDVI, но и идентификатор снимка 'imageID'. Последний позволяет не только различать снимки — из него можно извлечь дату съемки:

var imageID = image.id();
var date = ee.String(imageID).slice(0,8);

Есть и другое решение: можно составить мозаику для каждой даты съемки и проблема с несколькими значениями NDVI за одну дату исчезнет. Мозаики мы строить умеем (или здесь). Однако коллекция мозаик отнимает слишком много памяти, и потому позволит работать с меньшим числом участков, чем рассматриваемый нами подход.

#GEE

GEE-35. Временные ряды
Окончание. Начало

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

Участки с малым покрытием данными стоит исключить из рассмотрения (например, если данными NDVI покрыто меньше 10% площади участка), а в случае, когда участок имеет несколько NDVI за одну дату, доля покрытия данными может выступать в качестве весового коэффициента при расчете итогового среднего. Так или иначе, но мы должны получить долю покрытия данными каждого участка.

Площадь участка вычислим по полигону его границ:

var area = feature.area({'maxError': 1});

Число пикселей участка, покрытых данными, вычислим редьюсером ee.Reducer.count(). Функция combine позволяет объединить вычисление двух редьюсеров — mean и count:

...
reducer: ee.Reducer.mean().setOutputs(['ndvi']).combine({
      reducer2: ee.Reducer.count(),
      sharedInputs: true
    }), 
...    

Умножим число пикселей покрытых данными 'count', на площадь одного пикселя, и разделим на площадь участка area. Получим долю площади участка, покрытую данными:

var share = ee.Number(feature.get('count')).multiply(10).multiply(10).divide(area);

Кроме того, в результате нужно вернуть идентификатор участка 'id'. В итоге, возвращаем следующий набор свойств:

...
return feature.set({'id': feature.id(), 'ndvi': ndvi, 'imageID': imageID, 'date':date, 'share': share});

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

В результате мы получили временные ряды NDVI, которые нуждаются в доработке. Она зависит от того, как вы будете считать NDVI за одну дату, и какие снимки, в зависимости от доли покрытия данными, будете учитывать.

Можно ли вычислить готовые ряды NDVI полностью в GEE? Можно. Пример такого расчета приведен здесь. Но по нашему опыту подобное “допиливание” проще выполнять специализированным инструментом. Например, R или Python.

#GEE

Временной ряд, полученный в результате выполнения примера.

Космическая компания "Бюро 1440" выполнила сеанс видеосвязи из точки, где нет покрытия сотовой связи и земного интернета благодаря своей низкоорбитальной группировке.
Молодцы!

@dobriy_ovchinnikov

NRO ищет поставщиков новых видов спутниковых снимков

Национальное управление военно-космической разведки США (NRO) запрашивает предложения от коммерческих компаний, занимающихся дистанционным зондированием, на поставку новых перспективных видов данных оптического и инфракрасного диапазонов, кроме тех, что уже включены в контракты на поставку оптических снимков, которые управление заключило в прошлом году с компаниями Maxar Technologies, BlackSky и Planet Labs.

В заявке, опубликованной 14 августа, предлагаются такие технологии, как*: ночная съемка, съемка за пределами Земли, мультиспектральная съемка, видеосъемка и съемка в инфракрасном диапазоне.

Заказ данных осуществляется в рамках реализации программы Strategic Commercial Enhancements BAA**. В рамках этой программы NRO уже подписало соглашения с несколькими поставщиками радарных, радиочастотных и гиперспектральных данных.

Из указанных в заявке видов съемки самым непонятным является мультиспектральная. Вроде бы все перечисленные выше компании занимаются как раз мультиспектральной съемкой, особенно это касается Maxar.

Съемка за пределами Земли — это съемка космических объектов, в первую очередь, чужих спутников — новый сектор дистанционного зондирования, который надеются захватить компании вроде Maxar. После изменений в процедуре лицензирования, объявленных NOAA в начале августа, подобные данные могут предлагаться на коммерческой основе.


* Цитируем по SpaceNews: “the government needs imagery products across the electro-optical spectrum, including but not limited to panchromatic imagery, multispectral imagery, nighttime imaging, non-Earth imaging, video, and spectrally diverse Imagery including short-wave infrared, mid-wave infrared  and long-wave infrared.”

** BAA означает Broad Agency Announcement — инструмент, используемый правительственными агентствами для поддержки фундаментальных и прикладных исследований и той части разработок, которая не связана с созданием конкретной системы или закупкой оборудования.

#война

В России разрабатывают двигатель для работы космических аппаратов на сверхнизких околоземных орбитах

Ученые из Московского авиационного института (МАИ) и Московского государственного университета (МГУ) ведут работу над созданием двигательной установки для космических аппаратов, предназначенных для работы на сверхнизких околоземных орбитах высотой около 200 км. Двигатель будет создавать тягу, используя атмосферные газы — азот и кислород, которые еще сохраняются на данной высоте. Продолжительность работы такой установки, в отличие от традиционных, не ограничена запасом топлива. Ее применение позволит снизить размеры и стоимость спутников для таких высот, а также уменьшить затраты на их выведение в космос.

Сверхнизкие околоземные орбиты (Very Low Earth Orbit, VLEO) — это орбиты, лежащие ниже орбиты МКС (около 400 км). Использование таких орбит дает множество преимуществ. Для дистанционного зондирования это позволит получать снимки более высокого разрешения при тех же характеристиках съемочной аппаратуры, что и у современных спутников, находящихся на более высоких орбитах (400–700 км).

Проблема в том, что остатки атмосферы на сверхнизких орбитах создают значительное сопротивление движению спутника, быстро снижающее высоту его орбиты. Поддержать ее может двигательная установка, но современные двигатели требует слишком большого запаса топлива, что сводит на нет преимущества от использования сверхнизких орбит. Нужен новый, более экономичный двигатель. Эту задачу и решают специалисты из МАИ и МГУ. Они предполагают, что разработка может быть доведена до стадии опытного образца в течение 5–7 лет.

Отметим, что в России существует еще один проект двигателя для сверхнизких околоземных орбит, работа на которым ведется в настоящее время.

#VLEO

Диксон и его окрестности 3 марта 2023 года

Читатели спрашивают: можно ли полярной ночью сделать снимки в диапазонах видимого, ближнего инфракрасного и коротковолнового инфракрасного излучения?
Ответ однозначный: нельзя. Более того, снимков Sentinel-2 территорий, где наступает полярная ночь (ноябрь — февраль), нет в общем доступе. Скорее всего, в этот период съемка в этих районах не ведется.

А вот в период сумеречных ночей съемка возможна. Один из первых в этому году снимков порта Диксон сделан 3 марта. Все устье Енисейского залива покрыто сплошными льдами. Поселок находится в центре сцены. На острове Диксон темными линиями выделяется аэродром, значит снег на полосе убран. Немного южнее линии, соединяющей центры Диксона на острове и на материке, лежит зимняя дорога, соединяющая две части поселка. Ее можно разглядеть только на снимке без сжатия, прилагаемом ниже. Гораздо более заметный зимник лежит еще южнее, и тянется с востока на юг.

#севморпуть

GEE-36. Оценка облачного покрытия территории

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

Районом интереса являются окрестности порта Диксон 1️⃣. Оценивать долю облачного покрытия будем по среднемесячному значению свойства 'CLOUDY_PIXEL_PERCENTAGE' снимков Sentinel-2.

Алгоритм оценки состоит из трех шагов:

1. Задание района интереса.
2. Формирование коллекции снимков.
3. Выделение снимков заданного месяца и оценка среднего облачного покрытия для них.

Два первых шага нам хорошо знакомы. При формировании коллекции нужно взять данные за несколько лет, чтобы получить представительное количество снимков (мы брали за 3 года).

Будем выделять из коллекции снимки, относящиеся к определенному месяцу, и оценивать для них среднее значение свойства 'CLOUDY_PIXEL_PERCENTAGE'.

Для объединения (агрегации) снимков (Image) или векторных данных (Feature) в GEE существуют редьюсеры (Reducer), а для объединения свойств снимков и векторных данных используются агрегаторы — функции вида aggregate_что-нибудь, аргументом которых служит выбранное свойство. Например, сейчас нам понадобится ee.ImageCollection.aggregate_mean(), которая вычисляет среднее значение выбранного свойства среди объектов коллекции снимков. Вот как это выглядит:

var months = ee.List.sequence(1, 12);

var cloudList = months.map(function(month) {
  return col.filter(ee.Filter.calendarRange(month, month, 'month'))
            .aggregate_mean('CLOUDY_PIXEL_PERCENTAGE');
});

Результат показан на рисунке 2️⃣. Видно, что самыми малооблачными месяцами в данном районе являются февраль и март.

По опыту мы знаем, что получить хорошие снимки можно и в июле. Чтобы оценить количество малооблачных снимков за период, отберем из коллекции снимки, не превосходящие заданного порога облачного покрытия (например, 10%) и оценим число таких снимков за месяц, заменив агрегатор aggregate_mean на aggregate_count. Результат показан на рисунке 3️⃣. Таким образом, в июле за последние три года было получено 8 снимков с долей облачного покрытия менее 10%.

Существует около полутора десятков различных агрегаторов, вычисляющих как отдельные статистики, так и несколько статистик одновременно. Есть, в частности, aggregate_histogram, c помощью которого можно построить гистограмму снимков по доле облачного покрытия.

Код примера

#GEE