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

Эта статья описывает процесс подготовки к работе данных MOD14A1 (подробнее про эти данные). Исходно эти данные распространяются в HDF, который нельзя назвать удобным форматом для работы. Статья является синтезом других статей уже размещенных на нашем ресурсе, мы постараемся не дублировать информацию из них, а дать ссылки. Обратите внимание, что вам может быть не нужно осуществлять импорт именно в формат ESRI shape. Целевой формат зависит от того, что вам нужно сделать с данными в дальнейшем.

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

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

Для начала работы необходимо получить фрагменты (тайлы) MODIS на вашу территорию и временной промежуток. Данные в исходном формате HDF можно скачать здесь.

Для загрузки можно использовать скрипт (скачать). Строка запуска:

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

Начало работы

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

• Среда с установленным GDAL, python-gdal. Подойдет NextGIS QGIS.
• wget с дополнительными библиотеками
• Набор скриптов оберток для GDAL - DHI.

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

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

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

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

Конвертацию удобно проводить с помощью скрипта prepare_data.py (скачать). Строка запуска:

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

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

gdalinfo filename.hdf

Альтернативой prepare_data.py является MODIS Reprojection Tool (MRT). Его использование описано в статье: (Импорт продуктов MODIS уровней 2G, 3, 4 с помощью MRT).

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

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

Результатом предыдущего этапа является много растровых файлов в формате GeoTIFF. Если вам нужны именно растры, на этом можно остановиться. Если нет, конвертируем данные в вектор.

Схема работы следующая:

1. Конвертирует GeoTIFF в GRID
2. Конвертирует GRID в COVERAGE
3. Вытаскивает категории 8 и 9 (очаги пожаров с нормальной и повышенной достоверностью) в отдельные покрытия
4. Конвертирует эти покрытия в shape-файлы
5. Удаляет источники

Делать это удобно с помощью gdal_polygonize.py

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

Для пакетного склеивания ежедневных shape-файлов удобно использовать расширение Merge Shapes из QGIS (подробнее про расширение). Склейка в данном случае осуществляется потайлово, но можно склеивать и все вместе.

импорт

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

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

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

импорт

То есть используется 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 (источник) в ArcGIS (трансформация WGS84 -> Sinusoidal "на лету") приводит к достаточно хорошему соответствию, то есть система координат растра верна:

импорт

Для перепроецировании растра или вектора в 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

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