Зачастую для проведения операций растровой алгебры над несколькими растрами одновременно, необходимо, чтобы эти растры имели одинаковые экстент и разрешение. В данной статье описан процесс унификации экстента и разрешения нескольких растров при помощи скрипта, интегрируемого в Processing - модуль геопроцессинга QGIS.

[via Misanthrope's Thoughts]

[править] Установка скрипта

Загрузите архив со скриптом и help-файлом. Распакуйте архив и скопируйте содержимое папки "scripts" в директорию, предназначенную для Python-скриптов модуля Processing в QGIS (например, ~./qgis2/processing/scripts, если вы используете Linux). Если вы не знаете, где находится нужная папка перейдите в Processing -> Options and configuration -> Scripts -> Scripts folder, см. скриншот:

[править] Работа со скриптом

Перезапустите QGIS, если она была запущена. Откройте Processing Toolbox. Во вкладке Scripts вы увидите раздел Raster processing, а в нём - скрипт Unify extent and resolution:

Запустив его вы увидите диалоговое окно:

Обратите внимание, что во вкладке Help находится описание скрипта (на английском языке).

В поле "rasters" выберите 2 или более растра, экстент и разрешение которых надо сделать одинаковыми. Системы координат всех растров должны быть идентичными.

В поле "replace No Data value with" введите значение, которое будет использовано для новых пикселей, добавляемых в растр, а также заменит в результирующих растрах значения No Data обрабатываемых растров. Обратите внимание, что этому значению в результирующих растрах не будет присвоен атрибут No Data - они будут обрабатываться любым ПО, как обычные пиксели.

В поле "output directory" укажите путь к папке в которую должны быть сохранены результаты. Результирующие файлы будут названы следующим образом. Первая часть имени файла будет соответствовать названию исходного файла, к нему будет добавлен постфикс "_unified", например: "raster_1.tif" -> "raster_1_unified.tif".

После окончания основной работы скрипта появится сообщение о том, что растры были унифицированы и пользователю будет предоставлена возможность добавить результаты на карту, если он этого пожелает:

Для запуска скрипта в режиме пакетной обработки кликните по скрипту в Processing Toolbox правой кнопкой мыши и выберите "Execute as batch process":

Появится такое диалоговое окно, где вы сможете настроить все параметры, указанные выше:

[править] Как это работает

Опишем только основные идеи, использованные в алгоритме. Исходный код здесь подробно обсуждаться не будет, так как скрипт содержит свыше десятка функций и несколько классов (для Qt-сообщений). Кроме того, скрипт снабжен достаточным количеством комментариев.

Для обработки растров используется не QGIS API, а GDAL API. Работа алгоритма состоит в следующем:

• Для каждого угла обрабатываемых растров вычисляются значения X и Y. Из них выбираются максимальные и минимальные значения. Вычисление координат углов базируется на основе аффинной трансформации, лежащей в основе модели данных GDAL:

   Xgeo = GT[0] + Xpixel*GT[1] + Yline*GT[2]
   Ygeo = GT[3] + Xpixel*GT[4] + Yline*GT[5]

где Xgeo - значение координаты Х, Ygeo - значение координаты Y, Xpixel - номер столбца пикселя в растре, Ygeo - номер строки пикселя в растре, GT[2] и GT[4] вращение пикслея по осям X и Y соответственно, GT[1] - ширина пикселя, GT[5] - высота пикселя, GT[0] и GT[3] - соответственно X-координата и Y-координата верхнего левого угла верхнего левого пикселя ратсра, GT - параметры геотрансформации ратра, возвращаемые функцией gdal.GetGeoTransform().

• Далее берётся первый растр из поступивших на обработку и из него берётся информация о проекции и горизонтальных и вертикальных размерах пикселя. Эта информация вкупе с данными о финальном экстенте будет использована для создания болванок унифицированных растров.
• Наконец, для каждого из обрабатываемых растров создаётся унифицированная болванка на которую переносятся значения исходного растра, а оставшееся место (ведь экстент результирующего растра не меньше экстента исходного растра) забивается значениями, указанными пользователем. Это происходит следующим образом. Для каждого пикселя болванки рассчитываются её пространственные координаты (на основе формулы, представленной выше). Вычисляется, какому пикселю оригинального растра соответствуют данные координаты. Полученное значение пикселя оригинального растра записывается в болванку. Если координаты выходят за пределы оригинального растра или значение данного пикселя помечено, как No Data, то в болванку записывается значение, определённое пользователем.

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

Если вы таки захотите покопаться в исходном коде, то начните с функции main():

def main(rasters_list, no_data):
 
  # проверяем корректность указанного пути сохранения результатов
  if not checkFolder(output_directory):
    return NoPathGiven()
 
  # проверяем, совпадают ли проекции растров
  if not checkCRS(rasters_list):
    return WrongCRS()
 
  # получаем экстент для результирующих ратров
  fin_coordinates = finCoordinates(rasters_list)
 
  # собираем необходимую информацию для обработки
  processing_list = processingList(rasters_list, output_directory)
  raster_1 = rasters_list[0]
  raster_1 = gdal.Open(raster_1)
  main_geo_transform = raster_1.GetGeoTransform() # параметры геотрансформации эталонного растра
  proj = raster_1.GetProjection()
  if not no_data:
    no_data = raster_1.GetRasterBand(1).GetNoDataValue()
    if not no_data:
      no_data = -9999
 
  # обрабатываем растры
  for tupl in processing_list:
    raster = gdal.Open(tupl[0])
    output = ExtendRaster(raster, fin_coordinates, tupl[1], main_geo_transform, proj, no_data)
    raster = None
    if output != True:
      return None
 
  # выводим сообщение об успешном завершении
  result = Success()
  return result.initUI(processing_list)

[править] Вместо заключения

Пожелания и предложения по работе скрипта можно оставлять в соответствующей теме форума.

[править] Ссылки по теме

• Работа с растрами при помощи GDAL и Python
• Геопроцессинг с SEXTANTE для QGIS
• Подсчёт уникальных значений в растре при помощи SEXTANTE и QGIS