Практика использования raster tools

Эта страница является черновиком статьи.

В статье рассматриваются примеры практического использования набора растровых инструментов raster_tools.

Знакомство с raster_tools

В течении продолжительного времени развиваю инструментарий для препарирования георастров. Который является набором рецептов упрощающих взаимодействие с гибкой но сложной структурой python gdal. К сегодняшнему времени rater_tools напоминает некий «швейцарский нож» с разнообразным функционалом. Код проекта расположен на https://github.com/oldbay/raster_tools , самый простой способ его установки в каталоге пользователя - при помощи pip:

pip2.7 install --user git+https://github.com/oldbay/raster_tools

raster_tools состоит из следующих классов:

raster2array - Базовый класс, отвечающий за весь функционал загрузки растра и сохранение numpy массива выбранного канала(band). Выгруженный из растра массив может быть преобразован методами данного класса.

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

raster2transform(raster2array) - Субкласс, отвечающий за трансформацию растра: изменение размерности массива и перепроецирование. Так же есть поддержка методов базового класса.

raster2calc(raster2array) - Класс итерационного калькулятора. Он не наследует базовый класс raster2array, но использует некоторые его методы опосредованно.

raster2multiarray(raster2array) - Класс обработки мультирастров. Позволяет миксовать каналы мультиканального растра и обрабатывать его некоторыми методами raster2array.

multiarray2multiraster - Класс сохранения массивов, обработанных методами raster2multiarray, в мультиканальный растр.

Основной смысл существования всего этого «зоопарка» классов - это загрузка, обработка и выгрузка двухмерного numpy массива георастров. Помимо объекта 2-х мерного numpy массива некоторые методы raster2array позволяют выгружать данные в формате самоназвоного «стандартного словаря», имеющего следующий вид:

{
    "array": объект numpy.ndarray,
    "shape": кортеж (rows(Y), cols(X)) = numpy.ndarray.shape = (gdal.dataset.RasterYSize(), gdal.dataset.RasterXSize()),
    "transform": кортеж gdal.dataset.GetGeoTransform(),
    "projection": объект gdal.dataset.GetProjection()
}

Эта конструкция несёт в себе всю необходимую информацию для создания нового георастра, отличного от исходного.

Углубляться в подробное описание стрктуры классов и методов raster_tools не вижу необходимости - достаточно полная базовая документация есть в README.rst на https://github.com/oldbay/raster_tools . В данной статье хотелось бы показать «боевое» применение raster_tools на практике, в том виде как сам его в проектах «готовлю».

Для этих экспериментов был создан gis репозиторий с тестовыми данными , загрузить которые можно командой:

git clone  https://gitlab.com/oldbay/raster_tools_examples.git

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

git checkout origin/result

Тестовый репозирорий состоит из:

• Корень: содержит тестовые примеры - выполнение которых необходимо производить из текущего каталога;
• Каталог data: содержит исходные данные для экспериментов;
• Каталог logs: содержит результаты тестов расхода памяти и времени выполнения (branch result);
• Каталог module: содержит дополнительные python модули, необходимые для выполнение некоторых тестовых скриптов(в статье не описаны);
• Каталог results: содержит результаты экспериментов (branch result).

Растровый калькулятор.

Первоначально raster_tools был частью растрового калькулятора, сейчас отдельный проект raster_calc, поэтому исторически начнём с вычислений.

Пример простого растрового калькулятора представлен в calc.py (list 1) . В данном скрипте мультирастр разбирается на спектральные каналы, после чего вычисляется вегетативный индекс на основе rgb - TGI(украдено отсюда), Источником данных является data/multi.tif (img 1).

img 1 - data/multi.tif

list 1 - calc.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array, array2raster
import numpy as np

in_file = "data/multi.tif"
out_file = "result/calc.tif"

# загрузка каналов
red = raster2array(in_file, 1)
green = raster2array(in_file, 2)
blue = raster2array(in_file, 3)

# сохранение numpy массивов каналов в переменные.
r = red()
g = green()
b = blue()

# вычисление TGI
calc = np.choose(
    np.not_equal(g-r+b-255.0, 0.0),
    (
        -9999.0,
        np.subtract(
            g,
            np.multiply(0.39, r),
            np.multiply(0.61, b)
        )
    )
)

# сохранение вычисленного массива в растр.
array2raster(red, calc, out_file)

В результате выполнения calc.py создается растр result/calc.tif (img 2).

img 2 - result/calc.tif

У метода есть существенный минус - высокое потребление памяти. Фактически в памяти , на пике потребления, присутствуют сразу 4 массива растров: каналов r, g, b и вычисляемого индекса. При работе с большими георастрами это часто неприемлемо. Поэтому был разработан итерационный метод работы растрового калькулятора. Данный метод работает с растром посекторно в цикле. При таких вычислениях сильно снижается потребление памяти, но падает производительность.

Итерационный калькулятор реализован в calc_iter.py (list 2). В отличии от простого калькулятора здесь формула вычислений(в формате lambda) и переменные для вычислений(в виде объектов raster2array) передаются классу raster2calc.

list 2 - calc_iter.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array, array2raster, raster2calc
import numpy as np

in_file = "data/multi.tif"
out_file = "result/calc_iter.tif"

# загрузка каналов
red = raster2array(in_file, 1)
green = raster2array(in_file, 2)
blue = raster2array(in_file, 3)

# описание формулы вычисления TGI
calc_func = lambda r,g,b: np.choose(
    np.not_equal(g-r+b-255.0, 0.0),
    (
        -9999.0,
        np.subtract(
            g,
            np.multiply(0.39, r),
            np.multiply(0.61, b)
        )
    )
)

# вычисление TGI
# и сохраниение результата в формате
# стандартного словаря
calc = raster2calc()
out = calc(
        calc_func,
        r=red,
        g=green,
        b=blue
    )

# сохраниение стандартного словаря в растр.
array2raster(None, out, out_file)

raster2calc возвращает «стандартный словарь» преобразуемый классом array2raster в файл result/calc_iter.tif.

Если выполнить calc.py и calc_iter.py в обёртке memory_profiler - получаем графики потребления памяти:

Видно что итерационный метод выигрывает в потреблении памяти, но теряет в скорости расчётов. На текущем примере потери скорости не очень существенны - но при массовой обработке небольших георастров, бывает выгоднее пожертвовать памятью.

Вырезание области растра.

Иногда растровые вычисления необходимо производить только с некоторой выделенной областью георастра. Для «фигурного» вырезания областей, ограниченных координатами точек или bbox-ом полигонов, у raster2array есть ряд методов:

• cut_area(координаты точек) - метод вырезает область растра по координатам (не менее 2-х точек).
• cut_shp_file(имя shp файла) - метод вырезает растр по полигонам первого слоя shp файла.
• cut_ogr_geometry(wkt, geojson, gml, wkb) — метод вырезает растр по полигону загруженной геометрии.

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

Стандартный словарь применим и в растровом калькуляторе, для подобных вычислений необходимо вернуть numpy массив по ключу «array» list 1.1

list 1.1

# загрузка каналов в формате стандартного словаря
# выгрузка аналогична методам обрезки, но для
# всего исходного растра
red = raster2array(in_file, 1).get_std_dict()
green = raster2array(in_file, 2).get_std_dict()
blue = raster2array(in_file, 3).get_std_dict()

# сохранение numpy массивов каналов в переменные.
r = red()["array"]
g = green()["array"]
b = blue()["array"]

# вычисление TGI (как в list 1)

........


# сохранение вычисленного массива в стандартный словарь
out_std_dict = red
out_std_dict["array"] = calc

# сохранение стандартного словаря в растр.
array2raster(None, out_std_dict, out_file)

Теперь рассмотрим работу каждого метода обрезки растра:

Метод cut_area.

Позволяет обрезать растр по набору координат точек. Применяется когда необходимо обрезать область растра по массиву точек - для дальнейшей обработки. Пример использования приведён в cut_coords.py (list 3), в качестве исходных данных используется result/calc.tif.

list 3 - cut_coords.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array, array2raster

in_file = "result/calc.tif"
out_file = "result/cut_coords.tif"

# вырезание области растра по координатам
# и сохранение в формате стандартного словаря
std_dict = raster2array(in_file).cut_area(
    (296153.369,7137678.937),
    (296203.959,7137570.986),
    (296256.938,7137645.476)
)

# сохранение стандартного словаря в растр
array2raster(None, std_dict, out_file)

Результат выполнения метода сохраняется в out: result/cut_coords.tif (img 3.1 ,img 3.2)

img 3.1 - отображение области на исходном растре

img 3.2 - result/cut_coords.tif

Метод cut_shp_file.

Часто необходимо обрезать растр по геометрии слоя shp файла. Пример выполнения такой задачи приведён в cut_shp.py (list 4). Помимо исходного георастра result/calc.tif , скрипт использует шаблона обрезки - data/cut.shp.

list 4 - cut_shp.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array, array2raster

shp_file = "data/cut.shp"
in_file = "result/calc.tif"
out_file = "result/cut_shp.tif"

# вырезание области растра по полигону shp файла
# и сохранение в формате стандартного словаря
std_dict = raster2array(in_file).cut_shp_file(shp_file)

# сохранеие стандартного словаря в растр
array2raster(None, std_dict, out_file)

Результат сохраняется в георастр result/cut_shp.tif (img 4.1, img 4.2).

img 4.1 - геометрия из shp файла на фоне исходного растра

img 4.2 - result/cut_shp.tif

Обратите внимание что растр обрезан непосредственно по полигону, а не только по координатам bbox-а геометрии.

Метод cut_ogr_geometry.

Применяется если необходима обрезка георастра по геометрия в формате wkt, geojson, gml или wkb. Наиболее типичный способ использования - это вырезание растра на по wkt геометрии возвращённой из postgis.

Для выполнения примера следует провести подготовку(все примеры работы с БД приведены для linux):

• Установить postgresql и расширение postgis.
• В СУБД создать роль из под которой возможен логин в БД для загрузки дампа:

# su - postgres
$ psql
create role gis with login password 'gis';
create database dok_example owner gis;
ctrl-d

• Имя пользователя, пароль и имя базы можно использовать произвольные - но тогда необходимо будет изменить в list 5 переменные: dbname, dbuser, dbpass.
• Создать расширение postgis в БД.

$ psql -d doc_example
create extension postgis
ctrl-d

• Развернуть дамп из data/crowns.sql

psql -d doc_example < data/crowns.sql

• Проверить наличие python модуля psycopg2 - драйвера postgresql.(необходим для интерфейса запросов modules.db_interface)

После успешного завершения подготовительного этапа можно выполнять тестовый пример cut_wkt.py (list 5)

list 5 - cut_wkt.py

#! /usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array, array2raster
from modules import psql

# загрузка растра в объект raster2array
in_file = "result/calc.tif"
out_file = "result/cut_wkt.tif"

# установка параметров подключения в БД
dbhost = "localhost"
dbname = "dok_example"
dbuser = "gis"
dbpass = "gis"

# параметры поиска полигона в БД
geom_table = "crowns"
geom_id = 55775

# запрос в БД возвращающий полигон в формате WKT
SQL = """
select ST_AsText(wkb_geometry)
from {0}
where ogc_fid = {1}
""".format(
    geom_table,
    geom_id
)
_psql = psql(
    dbhost=dbhost,
    dbname=dbname,
    dbuser=dbuser,
    dbpass=dbpass
)
_psql.sql(SQL)
wkt_geom = _psql.fetchone()[0]
_psql.close()

# вырезание области растра по полигону WKT
# и сохранение в формате стандартного словаря
std_dict = raster2array(in_file).cut_ogr_geometry(wkt_geom)

# сохранеие стандартного словаря в растр
array2raster(None, std_dict, out_file)

В результате выполнения примера должен быть создан result/cut_wkt.tif (img 5.1, img 5.2).

img 5.1 - отображение геометрии из postgis на фоне исходного растра

img 5.2 - result/cut_wkt.tif

Трансформация.

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

Решить подобною проблему призван класс raster2transform. При помощи него можно изменить размерность всего массива, или только его части. Для работы тестового примера resize.py (list 6) требуется установка matplotlib и графической среды пользователя, исходным георастром является result/cut_shp.tif.

list 6 - resize.py

#! /usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array, raster2transform
from matplotlib import pyplot as plt

in_raster = "result/cut_shp.tif"
out_raster = "result/resize.tif"

# вывод массива в matplotlib
def plt_show(obj):
    plt.title("show")
    plt.imshow(obj.array(), cmap='gray')
    plt.show()

# загрузка растра в raster2array
raster = raster2array(in_raster)
print raster.rows, raster.cols

# кроме того возможна работа с:
#raster = in_raster #именем растрового файла
#raster = raster2array(in_raster).get_std_dict() #стандартным словарём

# изменение размерности массива (row, col) в цикле
for row, col in [(253, 210), (125, 105), (300, 300)]:
    raster = raster2transform(raster, row, col)
    print raster.rows, raster.cols
    plt_show(raster)

# сохранение растра c последней размерностью массива
raster.save(out_raster)

В результате работы скрипта 3 раза отобразится окно matplotlib.pyplot с изображением массива растра в следующих размерностях: (253,210)-img 6.1, (125,105)-img 6.2, (300,300)-img 6.3. Последний массив сохранится в растр result/resize.tif (img 6.3).

img 6.1 - разрешение (253x210)

img 6.1 - разрешение (125x105)

img 6.1 - разрешение (300x300)

Мультирастры.

Для вычисления индексов удобно работать с одноканальными растрами - представляя каждый канал как отдельно взятую сущность. Но иногда требуется создавать и работать с многоканальными конструкциями. Для это создана пара мультиканальных классов:

• raster2multiarray - Класс умеющий изменять очерёдность и количество каналов, при этом преобразующий многоканальный растр в стандартный словарь расширенного типа. От обычного он отличается тем что в поле «array» сохраняется не двухмерный (x, y), а трёхмерный numpy массив с индексом (bands, x, y). Для определения типа массива добавлено поле «multi_type» - принимающее значения: None - с массивом уже указанного индекса и «cv» - с индексом (x, y, bands). Массив типа multi_type = «cv» аналогичен возвращаемому функцией cv2.imread. Кроме того у raster2multiarray имеется набор методов обрезки растра аналогичных raster2array.

• multiarray2multiraster - Класс сохраняющий многоканальный стандартный словарь в мультиканальный растр.

Предлагаю два варианта использования данной пары классов, исходным растром для примеров послужит data/multi.tif:

Изменение последовательности каналов.

Для данного примера в rgb мультирастре переставим местами 1-й и 3-й каналы. Тестовый скрипт представлен в multi_rebands.py (list 7).

list 7 - multi_rebands.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2multiarray, multiarray2multiraster

img_in = 'data/multi.tif'
img_out = 'result/multi_rebands.tif'

# загрузка многоканального в объект raster2multiarray
# и изменение последовательности каналов 1-3, 2-2, 3-1
img = raster2multiarray(img_in, 3, 2, 1)

# выгрузка многоканального стандартого словаря
img = img.get_std_dict()

# cохранение многоканального стандартого словаря
# в мультирастр
multiarray2multiraster(img_out, img)

Результатом выполнения будет создание result/multi_rebands.tif (img 7).

img 7 - result/multi_rebands.tif

Обработка муьтирастров в opencv.

Теперь усложним задачу: удалим альфа канал мультирастра, передадим массив на обработку opencv и сохраним полученный результат. Для этого примера необходимо установить python opencv 2 или 3. Тестовый скрипт представлен в multi_opencv.py (list 8).

list 8 - multi_opencv.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

import numpy as np

from raster_tools import raster2multiarray, multiarray2multiraster
import cv2

img_in = 'data/multi.tif'
img_out = 'result/multi_opencv.tif'

# загрузка многоканального в объект raster2multiarray
# для обработки через cv2.bilateralFilter необходима загрузка только 1,2,3 канала
img = raster2multiarray(img_in, 1, 2, 3)

# установка типа возвращаемого многоканального массива
# в формат opencv
img.multi_type = "cv"

# установка типа данных возвращаемого массива как uint8
# необходимо для обработки данных в opencv
img.codage = np.uint8

# выгрузка многоканального стандартого словаря
std_dict = img.get_std_dict()

# работа в OpenCV c многоканальным массивом
std_dict["array"] = cv2.bilateralFilter(std_dict["array"],9,75,75)

# сохранение многоканального стандартного словаря
# в мультирастр
multiarray2multiraster(img_out, std_dict)

Результат будет сохранён в result/multi_opencv.tif (img 8).

img 8 - result/multi_opencv.tif

Работа с точками.

Не очевидный, но временами очень важный, метод препарирования георастров. Постановка задачи: запись значений точек растра с определёнными координатами в csv, БД или в виде поля в слой shp файла. Реализация последнего примера представлена в point_data_load_array.py (list 9). Исходными данными здесь выступают: георастр result/calc.tif и векторный файл со слоем точек data/tops.shp.

list 9 - point_data_load_array.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array
import ogr
import shutil

raster = "result/calc.tif"

# копируем тестовый shp файл с точками в каталог result
_filename = "tops"
_indir = "data"
_outdir = "result"
for _ex in ["dbf", "prj", "qpj", "shp", "shx"]:
    shutil.copyfile(
        "{0}/{1}.{2}".format(_indir, _filename, _ex),
        "{0}/{1}.{2}".format(_outdir, _filename, _ex),
    )

shp_file = "{0}/{1}.shp".format(_outdir, _filename)

# загрузка растра в объект raster2array
raster = raster2array(raster)

# сохранение numpy массива в объекте raster2array
raster.np_array_load()

# загрузка shp файла
_shp = ogr.Open(shp_file, update=1)
layer = _shp.GetLayerByIndex(0)

# добавление в shp файл нового поля
# для записи выгруженных значений растра
new_field = "POINT_DATA"
if layer.FindFieldIndex(new_field, 1) < 1:
    field = ogr.FieldDefn(new_field, ogr.OFTReal)
    layer.CreateField(field)

# обход циклом точек shp файла
for geometry in layer:

    # определение координат точки векторного файла
    x, _, y, _ = geometry.GetGeometryRef().GetEnvelope()

    # возврат значения данных точки растра
    # по координатам точки векторного файла
    raster_data = raster.get_pixel_value(float(x), float(y))

    # сохранение выгруженных из растра данных
    # в новое поле shp файла
    geometry.SetField(new_field, raster_data)
    layer.SetFeature(geometry)

В ходе выполнения скрипта исходный векторный файл data/tops.shp копируется в result/tops.shp, в который добавляется поле «POINT_DATA». В добавленном поле сохраняются значения точек растра по координатам геометрий векторного слоя. При помощи полученного векторного файла можно будет производить интерполяцию значений поля «POINT_DATA»(пример дальше по тексту).

В list 9 есть одна «хитрость»:

# сохранение numpy массива в объекте raster2array
raster.np_array_load()

Такая конструкция загружает numpy массив растра непосредственно в экземпляр класса rater2array - что уменьшает время выполнения операций, увеличивая потребление памяти. Для сравнения: в репозитории есть аналогичный скрипт point_data_dont_load_array.py где данная строка отсутствует. Сравнение графиков выполнения обоих скриптов в memory_profiler выглядит следующим образом:

"Ремонт" георастров.

Иногда процедурная генерация георастров некоторыми утилитами порождает «мутантов». При их отображении в gis программах проблема не заметна - так как геоданные будут правильно спозицианированны. Но при выгрузке из такого георастра numpy массива - оказывается что он зеркально отображён относительно оси x или y. Примером утилиты с такой странной процедурной генерации является gdal_grid(сообщение о проблеме на форуме).

Чтобы не быть голословными, выполним весь цикл процедур: интерполяция данных из предыдущего примера result/tops.shp при помощи gdal_grid, проверка полученного процедурного растра на «зеркальность», исправление его при помощи raster_tools. В качестве исходных данных будут использованы: result/tops.shp (источник данных для интерполяции) и result/calc.tif (для сравнения). Обработку произведём скриптом repair.py (list 11). Для работы repair.py необходим модуль modules.grid_utils и сама утилита gdal_rid.

list 11 - repair.py

#!/usr/bin/python2
# -*- coding: utf-8 -*-

from raster_tools import raster2array, array2raster
from modules import grid_utils

out_dir = "result"
shp_file = "result/tops.shp"
field_name = "POINT_DATA"
raster_file = "result/calc.tif"
grid_file = "result/{}.tif".format(field_name)
repair_file = "result/repair.tif"

# Интерполяция на основе shp файла
grid_utils(raster_file, shp_file, field_name, out_dir)

# загрузка результатов интерполяции в raster2array
grid = raster2array(grid_file)

# проверка "валидности" массива в растре
print grid.is_valid()

# возврат "отремонтированного" растра в формате стандартного словаря
grid = grid.repair()

# запись исправленного растра
array2raster(None, grid, repair_file)

В результате работы скрипта были созданы 2 растра: result/POINT_DATA.tif - результат работы gdal_grid и исправленный растр result/repair.tif. Сравним их gdalinfo лог:

$ gdalinfo result/POINT_DATA.tif
...
Pixel Size = (0.640086707766337,0.633172658674583)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  296008.694, 7137566.459) ( 40d46'52.75"E, 64d18'11.32"N)
Lower Left  (  296008.694, 7137767.175) ( 40d46'51.76"E, 64d18'17.78"N)
Upper Right (  296387.625, 7137566.459) ( 40d47'20.87"E, 64d18'12.13"N)
Lower Right (  296387.625, 7137767.175) ( 40d47'19.88"E, 64d18'18.60"N)
Center      (  296198.159, 7137666.817) ( 40d47' 6.31"E, 64d18'14.96"N)
Band 1 Block=592x1 Type=Float64, ColorInterp=Gray

$ gdalinfo result/repair.tif
...
Pixel Size = (0.640086707766337,-0.633172658674583)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  296008.694, 7137767.175) ( 40d46'51.76"E, 64d18'17.78"N)
Lower Left  (  296008.694, 7137566.459) ( 40d46'52.75"E, 64d18'11.32"N)
Upper Right (  296387.625, 7137767.175) ( 40d47'19.88"E, 64d18'18.60"N)
Lower Right (  296387.625, 7137566.459) ( 40d47'20.87"E, 64d18'12.13"N)
Center      (  296198.159, 7137666.817) ( 40d47' 6.31"E, 64d18'14.96"N)
Band 1 Block=256x256 Type=Float64, ColorInterp=Gray

и в качестве контрольной пробы:

$ gdalinfo result/calc.tif
...
Pixel Size = (0.064274599999999,-0.064274599999986)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  296007.704, 7137768.705) ( 40d46'51.68"E, 64d18'17.83"N)
Lower Left  (  296007.704, 7137564.891) ( 40d46'52.68"E, 64d18'11.27"N)
Upper Right (  296388.595, 7137768.705) ( 40d47'19.94"E, 64d18'18.65"N)
Lower Right (  296388.595, 7137564.891) ( 40d47'20.95"E, 64d18'12.08"N)
Center      (  296198.149, 7137666.798) ( 40d47' 6.31"E, 64d18'14.96"N)
Band 1 Block=256x256 Type=Float64, ColorInterp=Gray

Обратите внимание на второе значение Pixel Size(ось y): у растра сгенерированного gdal_grid это значение положительное, в отличии от исходного и исправленного георастров. Что и позволяет отображать POINT_DATA.tif в гис программах нормально, несмотря на наличие зеркально отображённого через ось y массива. То же самое сообщает нам вывод метода raster2array.is_valid():

# загрузка результатов интерполяции в raster2array
grid = raster2array(grid_file)
# проверка "валидности" массива в растре
print grid.is_valid()

Вывод строки:

[True, False]

Здесь ориентация растра по оси x истинное и ложное по оси y.

Послесловие.

Заканчивая описание способов употребления raster_tools хочу сказать что инструмент будет и дальше находится в состоянии постоянного развития. По мере поступления задач буду добавлять в него какие то новые приспособления, при этом постараюсь сохранить старые методы применения (чтобы не потерять совместимость с уже завершёнными проектами). Статья описывает состояние raster_tools версии 0.4 - на текущий момент последний stable. У raster_tools имеются следующие недостатки:

• проблемное перепроецирование векторных слоёв-шаблонов в проекцию растра, поэтому стараюсь использовать единую проекцию;
• обрезка по shp файлам «прибита гвоздями» к первому векторному слою;
• инструмент рассчитан на geotif - работа с другими растровыми форматами файлов в зачаточном состоянии;
• недостаточно оптимизирована работа с мультирастрами (не часто пока требуются).

В остальном инструмент меня устраивает, кого заинтересовал — отвечу на возникшие вопросы.