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

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

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

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

В течении продолжительного времени развиваю инструментарий для препарирования георастров. Который является набором рецептов, упрощающих взаимодействие с гибкой, но сложной структурой Python + GDAL. К сегодняшнему времени raster_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 на практике.

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

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 - получаем графики потребления памяти:

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

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

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

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

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

Стандартный словарь применим и в растровом калькуляторе, для подобных вычислений необходимо вернуть 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

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

Метод cut_ogr_geometry

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

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

• Установить PostgreSQL c расширением 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-массив растра непосредственно в экземпляр класса raster2array, что уменьшает время выполнения операций, увеличивая потребление памяти. Для сравнения: в репозитории есть аналогичный скрипт point_data_dont_load_array.py, где данная строка отсутствует. Сравнение графиков выполнения обоих скриптов в memory_profiler выглядит следующим образом:

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

Иногда процедурная генерация георастров некоторыми утилитами порождает «мутантов». При их отображении в ПО ГИС проблема не заметна - так как геоданные будут правильно позиционированы. Но при выгрузке из такого георастра 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_grid.

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 - на текущий момент последняя стабильная версия.

У raster_tools имеются следующие недостатки:

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