Распаковка информации о качестве данных MODIS: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
 
(не показано 39 промежуточных версий 2 участников)
Строка 1: Строка 1:
заготовка
{{Статья|Опубликована|modisqa}}
{{Аннотация|В статье разбирается вопрос распаковки и расшифровки информации о качестве данных продуктов MODIS}}


==Теория==
==Введение==
MODIS - камера дистанционного зондирования на борту спутников Terra и Aqua, снимающая каждую точку на Земле два раза в день и производящая десятки террабайт данных ежедневно. Данные MODIS проходят интенсивную обработку, включающую контроль качества. Помимо всего прочего, большое количество проблем доставляют облака, которые необходимо определенным образом маскировать. Поэтому все [http://gis-lab.info/qa/modislandprod.html продукты на базе данных MODIS] включают слой качества (он же QA=Quality Assessment, QC=Quality Control).
MODIS камера дистанционного зондирования на борту спутников Terra и Aqua, снимающая каждую точку на Земле два раза в день и производящая десятки терабайт данных ежедневно. Данные MODIS проходят интенсивную обработку, включающую контроль качества на всех уровнях продуктов. Помимо всего прочего, большое количество проблем, например, доставляют облака, которые необходимо определенным образом маскировать. Поэтому все [http://gis-lab.info/qa/modislandprod.html продукты на базе данных MODIS] включают слой качества (он же QA=Quality Assessment, QC=Quality Control), который рекомендуется использовать, если вы имеете дело с продуктами MODIS. Однако это не так просто из-за особой формы хранения данных.


Проблема: для экономии места данные QA "свернуты" в целочисленный вид, поэтому практически человеконечитаемы. Например, для продукта MCD13A2 значение <code>113</code> на самом деле значит следующее:  
Проблема: из целей практичного хранения создателями данные QA "свернуты" в целочисленный вид (UInt8, UInt16), поэтому практически человеконечитаемы. Чтобы их прочитать, нужно сначала их правильным образом распаковать. Например, для продукта MCD13A2 значение <code>113</code> на самом деле значит следующее:  


* Overall quality: Other Quality
* Overall quality: Other Quality
Строка 12: Строка 13:
* SCF_QC (five level confidence score): Main (RT) method failed due to problems other than geometry, empirical algorithm used
* SCF_QC (five level confidence score): Main (RT) method failed due to problems other than geometry, empirical algorithm used


Это статья посвящена тому, как распаковать данные MODIS QA в человекочитаемый вид для более обоснованного решения о фильтрации данных MODIS.
Эта статья посвящена тому, как распаковать данные MODIS QA в человекочитаемый вид для более обоснованного решения о выбраковке "плохих" пикселей. Сначала идет теоретическая часть рассказывающая, как устроены данные, а затем практика, показывающая, как сделать таблицу пересчета в LibreOffice Calc. По ряду продуктов приводятся готовые таблицы расшифровки QA.
 
==Теория==
Продукты MODIS распространяются в формате HDF, позволяющем хранить несколько растровых наборов данных (SUBDATASETS) одного разрешения. Кроме собственно данных, для которых сформирован этот продукт, в каждом HDF присутствует набор данных по оценке качества (QA — Quality Assessment).
 
[[Файл:Mcd15a3-h09v07-qa-legend.png|600px|thumb|center|Фрагмент растра QA продукта MCD15A3, h09v07 со значениями]]
 
===Набор данных QA===
Возьмем, например, один фрагмент данных по индексам растительности — [https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod13a2 MOD13A2] (Vegetation Indices 16-Day L3 Global 1km) и взглянем на него с помощью утилиты ''gdalinfo''.
 
<pre>gdalinfo MOD13A2.A2009209.h09v07.006.2015192134639.hdf</pre>
 
Эта команда выдает полное описание файла, включая все входящие в него наборы данных. Нас интересует последняя часть, там, где описаны поднаборы данных:
 
<syntaxhighlight lang="text">
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
</syntaxhighlight>
 
Более конкретно, нас интересует набор данных 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 относится к каждому пикселю.
 
Стоит обратить внимание, что наборов QA в продукте может быть несколько. См. например [https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd15a2 MCD15A2], где есть основной набор QA (FparLai_QC), а есть еще дополнительный (FparExtra_QC).
 
<syntaxhighlight lang="text">
  SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MCD15A2.A2006017.h00v08.005.2008063183615
.hdf":MOD_Grid_MOD15A2:FparLai_QC
  SUBDATASET_3_DESC=[1200x1200] FparLai_QC MOD_Grid_MOD15A2 (8-bit unsigned inte
ger)
  SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"MCD15A2.A2006017.h00v08.005.2008063183615
.hdf":MOD_Grid_MOD15A2:FparExtra_QC
  SUBDATASET_4_DESC=[1200x1200] FparExtra_QC MOD_Grid_MOD15A2 (8-bit unsigned in
teger)
</syntaxhighlight>
 
===Распаковка в биты===
 
При распаковке необходимо учитывать два основных момента.
 
1. Все продукты семейства HDF-EOS создаются в порядке '''от старшего к младшему''' ([https://ru.wikipedia.org/wiki/%D0%9F%D0%BE%D1%80%D1%8F%D0%B4%D0%BE%D0%BA_%D0%B1%D0%B0%D0%B9%D1%82%D0%BE%D0%B2 big-endian]). В этой схеме биты нумеруются справа налево.
 
Таким образом, первый бит в строке извлекается так:
 
<pre>=INT(MID(B13,8,1))</pre>
 
В таблицах, расшифровывающих QA ([https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod13a2 например]), это нулевой бит.
 
2. Дальше: каждый параметр QA может занимать 1 или несколько битов (битовое слово), и (важно!) нумерация битов в битовых словах идет наоборот — '''слева направо'''.
 
Например, число ''7425'' в двоичной системе исчисления записывается так (картинка [https://lpdaac.usgs.gov/sites/default/files/public/modis/docs/MODIS_LP_QA_Tutorial-1b.pdf отсюда], с исправлениями):
 
[[Файл:Modis-qa-bits-ru.png]]
 
===Расшифровка===
Для каждого продукта есть таблице, где описывается, какие биты QA что значат.
 
Расшифруем пример выше, имея в виду, что это продукт MOD09GQ, и его QA расшифрован в [https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod09gq этой таблице].
 
{| class="wikitable unsortable"
|-
! Битовое слово !! Описание
|-
| 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)
|-
|}


==Расчеты "на коленке"==
==Расчеты "на коленке"==
===План===
===План===
В качестве "коленки" будет использовать LibreOffice Calc. Вероятно те же самые формулы заработают и в MS Excel, но автор этой статьи это не проверял.
Наша задача — распаковать целочисленные коды из QA в человекочитаемый формат. Действовать будем по следующему плану:
 
Наша задача - распаковать целочисленные коды из QA в человекочитаемый формат. Действовать будем по следующему плану:


#Получим список уникальных значений из всех растров чтобы не пропустить чего-то полезного, но редкого.
#Получим список уникальных значений из всех растров чтобы не пропустить чего-то полезного, но редкого.
Строка 27: Строка 117:


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


Для каждого продукта есть страница описания, где можно подчерпнуть информацию о битах, их значениях и расшифровке. Вот например для [https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod13a2 MOD13A2].
Для каждого продукта есть страница описания, где можно почерпнуть информацию о битах, их значениях и расшифровке. Вот, например, для [https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod13a2 MOD13A2].


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


Например, для одного тайла MOD13A2:
Например, для одного тайла MOD13A2:
<pre>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</pre>
<pre>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</pre>


Вопросы извлечения большого количества данных, мозаицирования, перепроецирования в этой статье не рассматриваются.
Вопросы извлечения большого количества данных, мозаицирования, перепроецирования в этой статье не рассматриваются. Другие примеры использования инструментов GDAL в контексте работы с HDF можно посмотреть в этом [http://gis-lab.info/qa/gdal-examples.html#HDF4 наборе примеров].


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


Это легко делается с помощью скрипта на Python, который принимает на входе список папок, в которых надо обработать растры и имя текстового файла, куда будут сложены уникальные значений из них всех.
Это легко делается, например, с помощью [https://raw.githubusercontent.com/nextgis/dhi/master/general/get_unique.py скрипта на Python], который принимает на входе список папок или список растров, в которых надо обработать растры, и имя текстового файла, куда будут сложены уникальные значений из них всех.


<pre>python get_unique.py dir1,dir2,dir3 out.txt</pre>
<pre>python get_unique.py out.txt -fs dir1,dir2,dir3</pre>
или
<pre>python get_unique.py out.txt -rs raster1.tif,raster2.tif,raster3.tif</pre>


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


====Пересчет в биты====
Количество уникальных значений зависит от продукта, интенсивности его валидации, и может варьировать от десятков (MCD15A2) до сотен (MOD13A2).
Проблема в том, что по умолчанию фукнция DEC2BIN в Calc не работает с числами вне диапазона -512..512 (2^10 - один бит на четность). Поэтому, если битов качества больше чем 10, формулу придется модифицировать следующим образом, например для 16 bit:
 
====Пересчет в строку битов====
 
Пересчет в биты делается функцией BIN2DEC. Например, для 8 битного описания это делается так:
 
<pre>=DEC2BIN(A1,8)</pre>
 
Проблема в том, что по умолчанию функция DEC2BIN в Calc не работает с числами вне диапазона -512..512 (2^10 — один бит на четность). Поэтому, если битов в QA больше, чем 10, формулу придется модифицировать следующим образом, например, для 16 bit:


<pre>=DEC2BIN((MOD(A2,65536)/512),7)&DEC2BIN(MOD(A2,512),9)</pre>
<pre>=DEC2BIN((MOD(A2,65536)/512),7)&DEC2BIN(MOD(A2,512),9)</pre>
====Распаковка в слова и биты====
Создадим в таблице колонку для каждого бита и извлечем их из битовой строки простым MID-ом, учитывая порядок в строке и в битовом слове.
<pre>=INT(MID(B6,14,1))</pre>
===Расшифровка===
Можно нагородить IF-ов в формулах, но проще завести отдельный лист для словарей, где держать расшифровки кодов.
Например так, в первом столбце — битовое слово, во втором — краткая расшифровка (она попадет в таблицу), в третьем — полная:
<pre>VI usefulness
0000 Highest         Highest quality
0001 Lower         Lower quality
0010 Decreasing 1    Decreasing quality
0100 Decreasing 2    Decreasing quality
1000 Decreasing 3    Decreasing quality
1001 Decreasing 4    Decreasing quality
1010 Decreasing 5    Decreasing quality
1100 Lowest          Lowest quality
1101 Not useful      Quality so low that it is not useful
1110 L1B faulty      L1B data faulty
1111 Not processed  Not useful for any other reason/not processed
</pre>
Далее в столбцах расшифровки просто искать объединенные обратно в слова последовательности отдельных битов в этом словаре.
<pre>=VLOOKUP(CONCATENATE(E3,F3,G3,H3),Sheet2.$A$8:$C$18,2)</pre>
Словарей на отдельном листе таблицы может быть сколько угодно, только не забывайте менять координаты ячеек.


==Результаты==
==Результаты==
Здесь можно скачать готовые перекодировочные таблицы для следующих продуктов MODIS:
Готовые перекодировочные таблицы доступны в [https://github.com/nextgis/modis_quality_assessmet_tables репозитории Github].
 
На данный момент есть таблицы для следующих продуктов MODIS:


*MCD15A2 - Leaf Area Index - Fraction of Photosynthetically Active Radiation 8-Day L4 Global 1km
* MCD15A2: Leaf Area Index - Fraction of Photosynthetically Active Radiation 8-Day L4 Global 1km
*MOD13A2 - Vegetation Indices 16-Day L3 Global 1km
* MOD13A2: Vegetation Indices 16-Day L3 Global 1km. Эквивалент MOD13A1.
*MOD17A2 - Gross Primary Productivity 8-Day L4 Global 1km
* Gross Primary Productivity 8-Day L4 Global 1km
* MCD12Q1: Land Cover Type Yearly L3 Global 500 m SIN Grid.


[[Файл:Modis-qa-decode.png|750px|thumb|center|Фрагмент таблицы расшифровки QA продукта MOD13A2]]
[[Файл:Modis-qa-decode-vocab.png|690px|thumb|center|Фрагмент листа словарей QA продукта MOD13A2, используемых для расшифровки]]


==Полезные ссылки==
==Полезные ссылки==
* [https://lpdaac.usgs.gov/sites/default/files/public/modis/docs/MODIS_LP_QA_Tutorial-1b.pdf MODIS Land Products Quality Assurance Tutorial: Part 1]
* [https://lpdaac.usgs.gov/sites/default/files/public/modis/docs/MODIS_LP_QA_Tutorial-1b.pdf MODIS Land Products Quality Assurance Tutorial: Part 1]
* [https://lpdaac.usgs.gov/sites/default/files/public/modis/docs/MODIS_LP_QA_Tutorial-2.pdf MODIS Land Products Quality Assurance Tutorial: Part 2]
* [https://lpdaac.usgs.gov/sites/default/files/public/modis/docs/MODIS_LP_QA_Tutorial-2.pdf MODIS Land Products Quality Assurance Tutorial: Part 2]

Текущая версия от 15:17, 6 ноября 2016

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/modisqa.html


В статье разбирается вопрос распаковки и расшифровки информации о качестве данных продуктов MODIS

Введение

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

Проблема: из целей практичного хранения создателями данные QA "свернуты" в целочисленный вид (UInt8, UInt16), поэтому практически человеконечитаемы. Чтобы их прочитать, нужно сначала их правильным образом распаковать. Например, для продукта 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 в человекочитаемый вид для более обоснованного решения о выбраковке "плохих" пикселей. Сначала идет теоретическая часть рассказывающая, как устроены данные, а затем практика, показывающая, как сделать таблицу пересчета в LibreOffice Calc. По ряду продуктов приводятся готовые таблицы расшифровки QA.

Теория

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

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

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

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

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 относится к каждому пикселю.

Стоит обратить внимание, что наборов QA в продукте может быть несколько. См. например MCD15A2, где есть основной набор QA (FparLai_QC), а есть еще дополнительный (FparExtra_QC).

  SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MCD15A2.A2006017.h00v08.005.2008063183615
.hdf":MOD_Grid_MOD15A2:FparLai_QC
  SUBDATASET_3_DESC=[1200x1200] FparLai_QC MOD_Grid_MOD15A2 (8-bit unsigned inte
ger)
  SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"MCD15A2.A2006017.h00v08.005.2008063183615
.hdf":MOD_Grid_MOD15A2:FparExtra_QC
  SUBDATASET_4_DESC=[1200x1200] FparExtra_QC MOD_Grid_MOD15A2 (8-bit unsigned in
teger)

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

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

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

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

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

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

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

Например, число 7425 в двоичной системе исчисления записывается так (картинка отсюда, с исправлениями):

Modis-qa-bits-ru.png

Расшифровка

Для каждого продукта есть таблице, где описывается, какие биты 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 out.txt -fs dir1,dir2,dir3

или

python get_unique.py out.txt -rs raster1.tif,raster2.tif,raster3.tif

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

Количество уникальных значений зависит от продукта, интенсивности его валидации, и может варьировать от десятков (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)

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

Создадим в таблице колонку для каждого бита и извлечем их из битовой строки простым MID-ом, учитывая порядок в строке и в битовом слове.

=INT(MID(B6,14,1))

Расшифровка

Можно нагородить IF-ов в формулах, но проще завести отдельный лист для словарей, где держать расшифровки кодов.

Например так, в первом столбце — битовое слово, во втором — краткая расшифровка (она попадет в таблицу), в третьем — полная:

VI usefulness		
0000	Highest	        Highest quality
0001	Lower	        Lower quality
0010	Decreasing 1    Decreasing quality
0100	Decreasing 2    Decreasing quality
1000	Decreasing 3    Decreasing quality
1001	Decreasing 4    Decreasing quality
1010	Decreasing 5    Decreasing quality
1100	Lowest          Lowest quality
1101	Not useful      Quality so low that it is not useful
1110	L1B faulty      L1B data faulty
1111	Not processed   Not useful for any other reason/not processed

Далее в столбцах расшифровки просто искать объединенные обратно в слова последовательности отдельных битов в этом словаре.

=VLOOKUP(CONCATENATE(E3,F3,G3,H3),Sheet2.$A$8:$C$18,2)

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

Результаты

Готовые перекодировочные таблицы доступны в репозитории Github.

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

  • MCD15A2: Leaf Area Index - Fraction of Photosynthetically Active Radiation 8-Day L4 Global 1km
  • MOD13A2: Vegetation Indices 16-Day L3 Global 1km. Эквивалент MOD13A1.
  • Gross Primary Productivity 8-Day L4 Global 1km
  • MCD12Q1: Land Cover Type Yearly L3 Global 500 m SIN Grid.
Фрагмент таблицы расшифровки QA продукта MOD13A2
Фрагмент листа словарей QA продукта MOD13A2, используемых для расшифровки

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