GIS-Lab - Вклад [ru]

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T19:38:51Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
==== Импорт модулей ====
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
==== Описание аргументов ====
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
==== Счетчик дат ====
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 
==== Блок преобразования ====
Переменные id и od соответствуют исходной и конечной папкам, им присваиваются значения args.input_folder и args.output_folder. 
[[Файл:Ввод переменных.png]] 
Далее, для каждого канала растра из списка выполняются следующие операции: отбор нужных пикселей и их векторизация. В данном случае используется формат командной строки cmd. В первом случае мы используем утилиту gdal_calc.bat, где при помощи выражения «А*(А>7)» отбираем только те пиксели, значение которых строго больше 7, остальным присваиваем значение NoData. В результате из каждого исходного растра получаются 8 новых, названия которых соответствуют порядковому номеру дня в году. Затем, при помощи утилиты gdal_polygonize.bat преобразуем отобранные растровые элементы в формат shp-файла. Полученному файлу также присваивается порядковый номер дня в году, которому он соответствует. Команда os.system(cmd) запускает формат командной строки cmd. 
[[Файл:Блок_обработки_растра.png]] 
==== Блок присвоения атрибутивной информации ====
Следующий шаг в данном блоке позволяет присваивать атрибутивную информацию полученному shp-файлу. В начале создается поле и указывается его название и тип. Затем поле заполняется необходимой информацией. Например, в первом случае мы создаем текстовое поле «fulldate» с длинной до 9 символов и вносим в него год «str(b)» и порядковый номер дня в году «str(i)». Следующие два поля содержат в себе порядковый номер дня в году и год соответственно. 
[[Файл:Блок_заполнения_атрибутивной_таблицы.png]] 

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T19:36:42Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
==== Импорт модулей ====
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
Описание аргументов 
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
====Счетчик дат==== 
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 
====Блок преобразования==== 
Переменные id и od соответствуют исходной и конечной папкам, им присваиваются значения args.input_folder и args.output_folder. 
[[Файл:Ввод переменных.png]] 
Далее, для каждого канала растра из списка выполняются следующие операции: отбор нужных пикселей и их векторизация. В данном случае используется формат командной строки cmd. В первом случае мы используем утилиту gdal_calc.bat, где при помощи выражения «А*(А>7)» отбираем только те пиксели, значение которых строго больше 7, остальным присваиваем значение NoData. В результате из каждого исходного растра получаются 8 новых, названия которых соответствуют порядковому номеру дня в году. Затем, при помощи утилиты gdal_polygonize.bat преобразуем отобранные растровые элементы в формат shp-файла. Полученному файлу также присваивается порядковый номер дня в году, которому он соответствует. Команда os.system(cmd) запускает формат командной строки cmd. 
[[Файл:Блок_обработки_растра.png]] 
====Блок присвоения атрибутивной информации==== 
Следующий шаг в данном блоке позволяет присваивать атрибутивную информацию полученному shp-файлу. В начале создается поле и указывается его название и тип. Затем поле заполняется необходимой информацией. Например, в первом случае мы создаем текстовое поле «fulldate» с длинной до 9 символов и вносим в него год «str(b)» и порядковый номер дня в году «str(i)». Следующие два поля содержат в себе порядковый номер дня в году и год соответственно. 
[[Файл:Блок_заполнения_атрибутивной_таблицы.png]] 

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T19:34:18Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
=== Импорт модулей ===
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
Описание аргументов 
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
Счетчик дат 
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 
Блок преобразования 
Переменные id и od соответствуют исходной и конечной папкам, им присваиваются значения args.input_folder и args.output_folder. 
[[Файл:Ввод переменных.png]] 
Далее, для каждого канала растра из списка выполняются следующие операции: отбор нужных пикселей и их векторизация. В данном случае используется формат командной строки cmd. В первом случае мы используем утилиту gdal_calc.bat, где при помощи выражения «А*(А>7)» отбираем только те пиксели, значение которых строго больше 7, остальным присваиваем значение NoData. В результате из каждого исходного растра получаются 8 новых, названия которых соответствуют порядковому номеру дня в году. Затем, при помощи утилиты gdal_polygonize.bat преобразуем отобранные растровые элементы в формат shp-файла. Полученному файлу также присваивается порядковый номер дня в году, которому он соответствует. Команда os.system(cmd) запускает формат командной строки cmd. 
[[Файл:Блок_обработки_растра.png]] 
Блок присвоения атрибутивной информации 
Следующий шаг в данном блоке позволяет присваивать атрибутивную информацию полученному shp-файлу. В начале создается поле и указывается его название и тип. Затем поле заполняется необходимой информацией. Например, в первом случае мы создаем текстовое поле «fulldate» с длинной до 9 символов и вносим в него год «str(b)» и порядковый номер дня в году «str(i)». Следующие два поля содержат в себе порядковый номер дня в году и год соответственно. 
[[Файл:Блок_заполнения_атрибутивной_таблицы.png]] 

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T19:31:38Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
Импорт модулей 
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
Описание аргументов 
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
Счетчик дат 
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 
Блок преобразования 
Переменные id и od соответствуют исходной и конечной папкам, им присваиваются значения args.input_folder и args.output_folder. 
[[Файл:Ввод переменных.png]] 
Далее, для каждого канала растра из списка выполняются следующие операции: отбор нужных пикселей и их векторизация. В данном случае используется формат командной строки cmd. В первом случае мы используем утилиту gdal_calc.bat, где при помощи выражения «А*(А>7)» отбираем только те пиксели, значение которых строго больше 7, остальным присваиваем значение NoData. В результате из каждого исходного растра получаются 8 новых, названия которых соответствуют порядковому номеру дня в году. Затем, при помощи утилиты gdal_polygonize.bat преобразуем отобранные растровые элементы в формат shp-файла. Полученному файлу также присваивается порядковый номер дня в году, которому он соответствует. Команда os.system(cmd) запускает формат командной строки cmd. 
[[Файл:Блок_обработки_растра.png]] 
Блок присвоения атрибутивной информации 
Следующий шаг в данном блоке позволяет присваивать атрибутивную информацию полученному shp-файлу. В начале создается поле и указывается его название и тип. Затем поле заполняется необходимой информацией. Например, в первом случае мы создаем текстовое поле «fulldate» с длинной до 9 символов и вносим в него год «str(b)» и порядковый номер дня в году «str(i)». Следующие два поля содержат в себе порядковый номер дня в году и год соответственно. 
[[Файл:Блок_заполнения_атрибутивной_таблицы.png]] 

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T19:30:55Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
Импорт модулей 
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
Описание аргументов 
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
Счетчик дат 
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 
Блок преобразования 
Переменные id и od соответствуют исходной и конечной папкам, им присваиваются значения args.input_folder и args.output_folder. 
[[Файл:Ввод переменных.png]] 
Далее, для каждого канала растра из списка выполняются следующие операции: отбор нужных пикселей и их векторизация. В данном случае используется формат командной строки cmd. В первом случае мы используем утилиту gdal_calc.bat, где при помощи выражения «А*(А>7)» отбираем только те пиксели, значение которых строго больше 7, остальным присваиваем значение NoData. В результате из каждого исходного растра получаются 8 новых, названия которых соответствуют порядковому номеру дня в году. Затем, при помощи утилиты gdal_polygonize.bat преобразуем отобранные растровые элементы в формат shp-файла. Полученному файлу также присваивается порядковый номер дня в году, которому он соответствует. Команда os.system(cmd) запускает формат командной строки cmd. 
[[Файл:Блок_обработки_растра.png]] 
Блок присвоения атрибутивной информации
Следующий шаг в данном блоке позволяет присваивать атрибутивную информацию полученному shp-файлу. В начале создается поле и указывается его название и тип. Затем поле заполняется необходимой информацией. Например, в первом случае мы создаем текстовое поле «fulldate» с длинной до 9 символов и вносим в него год «str(b)» и порядковый номер дня в году «str(i)». Следующие два поля содержат в себе порядковый номер дня в году и год соответственно. 
[[Файл:Блок_заполнения_атрибутивной_таблицы.png]] 

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T19:29:13Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
Импорт модулей 
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
Описание аргументов 
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
Счетчик дат 
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 
Блок преобразования 
Переменные id и od соответствуют исходной и конечной папкам, им присваиваются значения args.input_folder и args.output_folder. 
[[Файл:Ввод переменных.png]] 
Далее, для каждого канала растра из списка выполняются следующие операции: отбор нужных пикселей и их векторизация. В данном случае используется формат командной строки cmd. В первом случае мы используем утилиту gdal_calc.bat, где при помощи выражения «А*(А>7)» отбираем только те пиксели, значение которых строго больше 7, остальным присваиваем значение NoData. В результате из каждого исходного растра получаются 8 новых, названия которых соответствуют порядковому номеру дня в году. Затем, при помощи утилиты gdal_polygonize.bat преобразуем отобранные растровые элементы в формат shp-файла. Полученному файлу также присваивается порядковый номер дня в году, которому он соответствует. Команда os.system(cmd) запускает формат командной строки cmd. 
[[Файл:Блок_обработки_растра.png]] 
Блок присвоения атрибутивной информации
Следующий шаг в данном блоке позволяет присваивать атрибутивную информацию полученному shp-файлу. В начале создается поле и указывается его название и тип. Затем поле заполняется необходимой информацией. Например, в первом случае мы создаем текстовое поле «fulldate» с длинной до 9 символов и вносим в него год «str(b)» и порядковый номер дня в году «str(i)». Следующие два поля содержат в себе порядковый номер дня в году и год соответственно. 
[[Файл:Блок_заполнения_атрибутивной_информации.png]] 

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Файл:Блок заполнения атрибутивной таблицы.png

2017-03-26T19:26:10Z

LeraBel:

Файл:Блок обработки растра.png

2017-03-26T19:25:07Z

LeraBel:

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T18:34:37Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
Импорт модулей 
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
Описание аргументов 
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
Счетчик дат 
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 
Блок преобразования 
Переменные id и od соответствуют исходной и конечной папкам, им присваиваются значения args.input_folder и args.output_folder. 
[[Файл:Ввод переменных.png]] 
Далее, для каждого канала растра из списка выполняются следующие операции: отбор нужных пикселей и их векторизация. В данном случае используется формат командной строки cmd. В первом случае мы используем утилиту gdal_calc.bat, где при помощи выражения «А*(А>7)» отбираем только те пиксели, значение которых строго больше 7, остальным присваиваем значение NoData. В результате из каждого исходного растра получаются 8 новых, названия которых соответствуют порядковому номеру дня в году. Затем, при помощи утилиты gdal_polygonize.bat преобразуем отобранные растровые элементы в формат shp-файла. Полученному файлу также присваивается порядковый номер дня в году, которому он соответствует. Команда os.system(cmd) запускает формат командной строки cmd.

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T18:31:51Z

LeraBel: /* Конвертация в формат ESRI shape */

Эта статья описывает процесс подготовки к работе данных MOD14A1 ([http://gis-lab.info/qa/mod14a1.html подробнее] про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

Многие этапы можно реализовать по другому. Исправления к процессу аналогичные по скорости и удобству всячески приветствуются.

==Получение данных==

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать [http://gis-lab.info/qa/mod14a1.html здесь].

Для загрузки можно использовать скрипт ([https://github.com/nextgis/dhi/blob/master/general/download_modis_data.py скачать]). Строка запуска:

<pre>python download_modis_data.py 2016 46 http://e4ftl01.cr.usgs.gov/MOLT/MOD14A1.006/ hdf</pre>

==Начало работы==
Для начала работы понадобится:

* Среда с установленным GDAL, python-gdal. Подойдет [http://nextgis.ru/nextgis-qgis NextGIS QGIS].
* wget с дополнительными библиотеками
* Набор скриптов оберток для GDAL - [https://github.com/nextgis/dhi DHI]. В принципе они не обязательны, но могут оказаться полезными для автоматизации.

==Конвертация в GeoTIFF==

Исходная проекция данных - Синусоидальная (Sinusoidal), конвертация в неё, а не в другую, возможно более удобную сэкономит время на перепроецирование. Перепроецировать в другую систему координат лучше окончательный векторный слой.

Исходное разрешение: 1000 метров.

Исходные наборы данных: Формат HDF позволяет хранить в одном файле несколько наборов данных (SDS, subdatasets). Нам нужно импортировать только данные за 8 дней (допустим каналы с оценкой качества - Quality Assessment пока не нужны).

Конвертацию удобно проводить с помощью скрипта prepare_data.py ([https://github.com/nextgis/dhi/blob/master/general/prepare_data.py скачать]). Строка запуска:

<pre>python prepare_data.py MODIS_Grid_Daily_Fire:FireMask с:\MCD14A1\2003\hdf\ с:\MCD15A1\2003\tif\ -e SIN</pre>

Узнать строку названия SDS (в примере выше это MODIS_Grid_Daily_Fire:FireMask) можно посмотрев метаданные одного из файлов HDF с помощью gdalinfo:

<pre>gdalinfo filename.hdf</pre>

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: ([http://gis-lab.info/qa/modisimport.html Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT]).

Количество отдельных растров генерируемых при импорте данных с 2000 по 2010 год составляет около 3950 растров.

==Конвертация в формат ESRI shape==

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.
Полученные файлы представляют собой 8-канальные растры (каждый канал соответствует одному дню), в которых каждой ячейке присвоено значение от 0 до 9.
{| class="wikitable"
|-
! Код ячейки !! Значение ячейки !! Значение ячейки
|-
| 0 || отсутствует исходная информация || missing input data
|-
| 1 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 2 || не обрабатывалось (устарело) || not processed (obsolete)
|-
| 3 || вода || water
|-
| 4 || облако || cloud
|-
| 5 || не пожар || non-fire
|-
| 6 || неизвестно || unknown
|-
| 7 || пожар (низкая достоверность) || fire (low confidence)
|-
| 8 || пожар (номинальная достоверность) || fire (nominal confidence)
|-
| 9 || пожар (высокая достоверность) || fire (high confidence)
|}
Далее необходимо получить данные о пожарах для каждого дня в году, путем разложения растра поканально и извлечением только тех ячеек, значение которых равно 8 или 9. Для удобства оперирования данными в дальнейшем, их необходимо преобразовать в формат shp-файла.Для разложения растры на каналы используется команда gdal_translate. Ниже приведен пример запуска: 
<pre>gdal_translate -b 1 input.tif output.tif</pre>
где -b указывает на то, что нужно извлечь канал, 1 - номер канала 
С помощью команды gdal_calc необходимо извлечь только те ячейки, значение которых равно 8 ли 9, что соответствует пожарам с номинальной и высокой степенью достоверности, а остальным значениям присвоить значение NoData. Условием отбора является выражение “А*(А>7)”. Пример запуска команды: 
<pre>gdal_calc -A input.tif --outputfile= output.tif --calc='A*(A>7)' --NoDataValue=0</pre>
С помощью команды gdal_polygonize можно преобразовать полученные ячейки в полигоны и присвоить им атрибутивную информацию в виде даты возгорания.
<pre>gdal_polygonize -nomask output.tif output.shp</pre>
где при необходимости можно указать файл маски, при этом поменяв параметр "-nomask" на "-mask" и искомый файл. 
Таким образом можно легко извлечь необходимую информацию и преобразовать её в формат ESRI shapefile.

Также можно воспользоваться скриптом proc_fields.py, позволяющим создавать shp-файлы, содержащие полигоны зафиксированных в течение 1 дня горячих точек из исходного 8-канального растра, полученного на предыдущем этапе. 
Для каждого 8-канального растра скрипт выполняет следующие действия: 
- внутри каждого канала исходного изображения отбирает только те ячейки, значение которых равно 8 (номинальная достоверность наличия пожара в данной точке) и 9 (высокая достоверность наличия пожара), всем остальным ячейкам присваивает значение 0 (NoData); 
- векторизует полученные данные с помощью команды gdal_polygonize. Результат: shp-файлы с полигонами пожаров, где название файла соответствует порядковому номеру дня в году; 
- добавляет в каждый сформированный shp-файл поле со значением соответствующей даты возгорания в виде порядкового номера дня в году (н-р, 1 января - "1", 31 декабря - "365"). 
В дальнейшем, при необходимости, можно объединить ряд файлов в один, тем самым сформировав файл пожаров за определенный год, месяц или неделю (см. пункт 4). 
Схема преобразования данных представлена на рисунке ниже: 
[[Файл:Схема скрипта.jpg]]
Рассмотрим подробнее структуру скрипта: 
Импорт модулей 
Для данного скрипта необходимы следующие модули: os, argparse и datetime из которого нам нужны классы date и timedelta. Класс datetime.date позволяет выделить из полного формата даты только год, месяц и день. Класс datetime.timedelta позволяет высчитать разницу между двумя моментами времени, с точностью до микросекунд. Модуль os включает множество функций для работы с операционной системой, причём их поведение, не зависит от ОС, поэтому программы остаются переносимыми. Модуль argparse позволяет обрабатывать опции и аргументы командной строки, которой выполняется скрипт. В данном случае можно обойтись и без него, но он позволяет вводить нужную информацию (год, пути исходного и конечного файлов) непосредственно при запуске из командной строки, тем самым исключая необходимость внесения правок в само тело скрипта. 
[[Файл:Пример_импорта_модулей.png|Пример импорта модулей]] 
Описание аргументов 
Следующим этапом является описание аргументов, которые, как было сказано выше, вводятся в командной строке. К ним относятся год «year», исходная папка с данными «input_folder» и папка «output_folder», в которую будут помещены результаты. 
[[Файл:Пример_описания_аргументов.png]] 
Строка «args = parser.parse_args()» является активатором запуска модуля parse. Далее, вводится переменная b, которой присваивается значение вводимого года исходных данных. 
Счетчик дат 
Исходные данные представляют собой 8-ми канальные растры, где каждый канал соответствует одному дню. Имя файла содержит в себе дату первого дня из восьмидневного периода. Например, «2016.01.01.tif», «2016.01.09.tif» и т.д. Таким образом, для того, чтобы скрипт работал, необходимо ввести имена всех исходных файлов в виде списка. Однако, данную процедуру необходимо будет повторять для каждого года в отдельности, особенно для високосных годов, где помимо года в названии файла меняются месяц и день. Счетчик дат необходим для того, чтобы не вводить вручную имена исходных растров, а автоматически создавать список этих имен. 
Необходимо создать пустой список, в который будут помещаться имена исходных файлов после их создания. В данном случае ему присваивается имя rasters_list. Введем переменную «а» и присвоим ей значение равное 0. При обработке каждого последующего из исходных 46 растров, переменная «a» будет увеличиваться на 8. Переменная delta отвечает за счет дат, т.е. прибавляет к каждой предыдущей дате переменную «а», соответствующую количеству дней, прошедших с момента отсчета (1 января каждого года) до начала очередного 8-ми дневного периода. Исходная дата «base_date» соответствует первому января каждого года. Соответственно, для того, чтобы получить следующую дату, необходимо прибавить к ней переменную delta. Следующим шагом является преобразование формата даты в строку, для этого используется класс time.strftime. Дате, представленной в формате строки, прибавляется окончание «.tif», соответствующее расширению растра, и список пополняется при помощи команды «append». 
[[Файл:Создание_и_запосление_списка_дат.png]] 

==4. Создание единого слоя==

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS ([http://gis-lab.info/qa/merge-shapes.html подробнее про расширение]). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

<center>[[Image:mod14a1-import-03.gif|393px|импорт]]</center>

==5. Перепроецирование в WGS84 Lat/Long==

Результирующие shape-файлы может понадобится перепроецировать из синусоидальной проекции в географическую систему координат WGS84. Можно делать это и на этапе импорта в GeoTIFF, но это не рекомендуется, так как займет во много раз больше времени чем перепроецировани уже отфильтрованного вектора.

Система координат растров после работы MRT следующая:

<center>[[Image:mod14a1-import-04.gif|677px|импорт]]</center>

То есть используется Sinusoidal проекция на сфере радиусом 6370997 метров (цифра подтверждается MRT_Users_Manual.doc, стр. 57). В ArcGIS такая сфера соответствует сфере D_Sphere_ARC_INFO из папки Geographic Coordinate Systems\Spheroid based\. Соответственно для векторного файла PRJ файл будет выглядеть следующим образом:

PROJCS["World_Sinusoidal",
GEOGCS["GCS_Sphere_ARC_INFO",
DATUM["D_Sphere_ARC_INFO",
SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]
],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]
],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]
]

Наложение растра только что импортированного MRT и границ субъектов РФ в Lat/Long WGS84 ([http://gis-lab.info/qa/rusbounds-rosreestr.html источник]) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

<center>[[Image:mod14a1-import-05.gif|638px|импорт]]</center>

Для перепроецировании растра или вектора в Lat/Long WGS84 используем:

gdalwarp -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 A2000057_2.tif A2000057_2_wgs.tif

или

ogr2ogr -s_srs "+proj=sinu +R=6370997.0 +nadgrids=@null +wktext" -t_srs EPSG:4326 h19v02_wgs.shp h19v02.shp

Стоит повторить, что это лишь один из возможных путей импорта данных в векторных формат, однако он работает.

Файл:Создание и запосление списка дат.png

2017-03-26T18:28:57Z

LeraBel:

Файл:Ввод переменных.png

2017-03-26T18:28:44Z

LeraBel:

Файл:Пример описания аргументов.png

2017-03-26T18:28:32Z

LeraBel:

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T18:26:30Z

LeraBel: /* Конвертация в формат ESRI shape */

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T18:25:40Z

LeraBel: /* Конвертация в формат ESRI shape */

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T18:24:00Z

LeraBel: /* Конвертация в формат ESRI shape */

Файл:Пример импорта модулей.png

2017-03-26T18:23:17Z

LeraBel:

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T18:17:16Z

LeraBel: /* Конвертация в формат ESRI shape */

Импорт данных MOD14A1 в формат ESRI shape

2017-03-26T18:15:32Z

LeraBel: /* Конвертация в формат ESRI shape */

Файл:Схема скрипта.jpg

2017-03-26T18:14:48Z

LeraBel:

Импорт данных MOD14A1 в формат ESRI shape

2016-10-06T17:13:34Z

LeraBel: /* Конвертация в формат ESRI shape */

Импорт данных MOD14A1 в формат ESRI shape

2016-09-22T13:21:17Z

LeraBel: /* Конвертация в формат ESRI shape */

Файл:Gdal translate.jpg

2016-09-22T13:07:51Z

LeraBel:

Импорт данных MOD14A1 в формат ESRI shape

2016-06-07T12:46:03Z

LeraBel: /* Конвертация в GeoTIFF */