Распаковка информации о качестве данных MODIS

заготовка

Введение

MODIS - камера дистанционного зондирования на борту спутников Terra и Aqua, снимающая каждую точку на Земле два раза в день и производящая десятки террабайт данных ежедневно. Данные MODIS проходят интенсивную обработку, включающую контроль качества. Помимо всего прочего, большое количество проблем доставляют облака, которые необходимо определенным образом маскировать. Поэтому все продукты на базе данных MODIS включают слой качества (он же QA=Quality Assessment, QC=Quality Control).

Проблема: для экономии места данные QA "свернуты" в целочисленный вид, поэтому практически человеконечитаемы. Например, для продукта MCD13A2 значение 113 на самом деле значит следующее:

• Overall quality: Other Quality
• Source: Terra
• Detectors: Detectors are OK
• CloudState: Mixed clouds
• SCF_QC (five level confidence score): Main (RT) method failed due to problems other than geometry, empirical algorithm used

Это статья посвящена тому, как распаковать данные MODIS QA в человекочитаемый вид для более обоснованного решения о фильтрации данных MODIS. Сначала идет теоретическая часть рассказывающая как устроены данные, а затем практика, показывающая как сделать таблицу пересчета в LibreOffice Calc.

Теория

Продукты MODIS распространяются в формате HDF, позволяющем хранить несколько растровых наборов данных (SUBDATASETS) одного разрешения. Кроме собственно данных для которых сформирован этот продукт, в каждом HDF присутствует набор данных по оценке качества (QA - Quality Assessment).

Возьмем например один фрагмент данных по индексам растительности - MOD13A2 (Vegetation Indices 16-Day L3 Global 1km) и взглянем на него с помощью gdalinfo.

Набор данных QA

gdalinfo MOD13A2.A2009209.h09v07.006.2015192134639.hdf

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

Subdatasets:
  SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD13A2.A2009209.h09v07.006.2015192134639
.hdf":MODIS_Grid_16DAY_1km_VI:1 km 16 days NDVI
  SUBDATASET_1_DESC=[1200x1200] 1 km 16 days NDVI MODIS_Grid_16DAY_1km_VI (16-bi
t integer)
  SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"MOD13A2.A2009209.h09v07.006.2015192134639
.hdf":MODIS_Grid_16DAY_1km_VI:1 km 16 days EVI
  SUBDATASET_2_DESC=[1200x1200] 1 km 16 days EVI MODIS_Grid_16DAY_1km_VI (16-bit
 integer)
  SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MOD13A2.A2009209.h09v07.006.2015192134639
.hdf":MODIS_Grid_16DAY_1km_VI:1 km 16 days VI Quality
  SUBDATASET_3_DESC=[1200x1200] 1 km 16 days VI Quality MODIS_Grid_16DAY_1km_VI
(16-bit unsigned integer)
  SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"MOD13A2.A2009209.h09v07.006.2015192134639
.hdf":MODIS_Grid_16DAY_1km_VI:1 km 16 days red reflectance

Более конкретно, нас интересует набор данных SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MOD13A2.A2009209.h09v07.006.2015192134639 .hdf":MODIS_Grid_16DAY_1km_VI:1 km 16 days VI Quality. Все что содержит в описании слово Quality является частью QA.

Из описания набора можно видеть, что он такого же размера как и основные растры (1 и 2 наборы ) - 1200х1200, из чего становится ясно, что QA относится к каждому пикселю.

Распаковка в биты

При распаковке необходимо учитывать два основных момента.

1. Все продукты семейства HDF-EOS создаются в порядке от старшего к младшему (big-endian). В этой схеме биты нумеруются справа-налево.

Таким образом, первый бит в строкеизвлекается так:

=INT(MID(B13,8,1))

В таблицах расшифровывающих QA (например) это нулевой бит.

2. Дальше, каждый параметр QA может занимать 1 или несколько битов (битовое слово) и, важно!, нумерация битов в битовых словах идет наоборот слева-направо.

Расшифровка

Для каждого продукта есть таблице, где описывается какие биты QA что значат.

Расшифруем пример выше имея в виду что это продукт MOD09GQ и его QA расшифрован в этой таблице.

Битовое слово	Описание
01	Less than ideal quality some or all bands
00	Clear (Cloud state)
0000	Highest quality (Band 1 data quality
1101	Correction out of bounds pixel constrained to extreme allowable value
1	Yes (Atmospheric correction performed)
0	No (Adjacency correction performed)
00	Spare (unused)

Расчеты "на коленке"

План

Наша задача - распаковать целочисленные коды из QA в человекочитаемый формат. Действовать будем по следующему плану:

1. Получим список уникальных значений из всех растров чтобы не пропустить чего-то полезного, но редкого.
2. Каждый код переведем в биты
3. В соответствии со схемой QA для конкретного продукта нарезаем всю последовательность на битовые слова
4. Раскладываем каждое слово на биты слева направо
5. Складываем обратно и используем перекодировочные таблицы для того чтобы "расшифровать" каждую последовательность.

Реализация

В качестве "коленки" будет использовать LibreOffice Calc. Вероятно те же самые формулы заработают и в MS Excel, но автор этой статьи это не проверял. Реализовывать распаковку будем в Calc с целью визуализировать и прочитать, осмыслить все значения. В будущем разумеется было бы неплохо придумать для этого скрипт на Python или библиотеку. Так как эта процедура уже была произведена для ряда продуктов (см. ниже), используем таблицу для одного из них в качестве шаблона.

Для каждого продукта есть страница описания, где можно подчерпнуть информацию о битах, их значениях и расшифровке. Вот например для MOD13A2.

Извлечение QA растров

Раньше для получения данных MODIS из HDF было принято использовать MRT, но сейчас извлечение растров QA стало удобно делать с помощью gdal_translate.

Например, для одного тайла MOD13A2:

gdal_translate HDF4_EOS:EOS_GRID:"MOD13A2.A2005001.h00v08.006.2015157042257.hdf":"MODIS_Grid_16DAY_1km_VI:1 km 16 days VI Quality" output.tif

Вопросы извлечения большого количества данных, мозаицирования, перепроецирования в этой статье не рассматриваются. Другие примеры использования инструментов GDAL в контексте работы с HDF можно посмотреть в этом наборе примеров.

Получение уникальных значений

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

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

python get_unique.py dir1,dir2,dir3 out.txt

Результат это просто столбец целочисленных значений, которые вставляются в нашу таблицу - шаблон.

Количество уникальных значений зависит от продукта, интенсивности его валидации и может варьировать от десятков (MCD15A2) до сотен (MOD13A2).

Пересчет в строку битов

Пересчет в биты делается функцией BIN2DEC. Например, для 8 битного описания это делается так:

=DEC2BIN(A1,8)

Проблема в том, что по умолчанию фукнция DEC2BIN в Calc не работает с числами вне диапазона -512..512 (2^10 - один бит на четность). Поэтому, если битов в QA больше чем 10, формулу придется модифицировать следующим образом, например для 16 bit:

=DEC2BIN((MOD(A2,65536)/512),7)&DEC2BIN(MOD(A2,512),9)

Результаты

Здесь можно скачать готовые перекодировочные таблицы для следующих продуктов MODIS:

• MCD15A2 - Leaf Area Index - Fraction of Photosynthetically Active Radiation 8-Day L4 Global 1km
• MOD13A2 - Vegetation Indices 16-Day L3 Global 1km
• MOD17A2 - Gross Primary Productivity 8-Day L4 Global 1km

Фрагмент растра QA продукта MCD15A3, h09v07 со значениями

Полезные ссылки

• MODIS Land Products Quality Assurance Tutorial: Part 1
• MODIS Land Products Quality Assurance Tutorial: Part 2