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

Материал из GIS-Lab
Перейти к навигации Перейти к поиску

Эта статья описывает процесс подготовки к работе данных 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

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

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

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

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

Так как это требуется нашему конвертеру в вектор, сначала переименовываем все растры так, чтобы имена стали покороче, например по шаблону: [N9-16]_[N72]. Таким образом из названия MOD14A1.A2000065.H19V01.005.2006270190019_res.FireMask.Number_of_Days_04.tif, мы получим более лаконичное: A2000065_b4.tif

Примечание: иногда в исходных файлах HDF содержится меньше чем 8 каналов, в этом случае будут импортированы дополнительные каналы не являющиеся пожарными масками (QA, MaxFRP, sample). После окончания процесса конвертации их рекомендуется удалить.

После этого запускаем скрипт для Arcinfo Workstation. Этот скрипт делает следующее:

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

Подробнее про запуск, содержание и работу с подобными скриптами можно ознакомиться в статье "Пакетный импорт растровых данных в GRID и их обработка с помощью Arcinfo Workstation" ([importraster-ai.html прочитать]).

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

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

импорт

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 ([rusbounds-rosreestr.html источник]) в 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

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

Источник — https://wiki.gis-lab.info/index.php?title=Импорт_данных_MOD14A1_в_формат_ESRI_shape&oldid=24166