Практика использования 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 - работа с другими растровыми форматами файлов в зачаточном состоянии;
• недостаточно оптимизирована работа с мультирастрами (не часто пока требуются).

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