Практика использования raster tools: различия между версиями
Oldbay (обсуждение | вклад) |
Oldbay (обсуждение | вклад) |
||
Строка 520: | Строка 520: | ||
</syntaxhighlight> | </syntaxhighlight> | ||
Результат будет сохранён в result/multi_opencv.tif. | Результат будет сохранён в result/multi_opencv.tif (img 8). | ||
img 8 - result/multi_opencv.tif | |||
[[Файл:Multi_opencv.png|500px]] | |||
==Работа с точками.== | ==Работа с точками.== |
Версия от 13:13, 1 июня 2018
В статье рассматриваются примеры практического использования набора растровых инструментов 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 - работа с другими растровыми форматами файлов в зачаточном состоянии;
- недостаточно оптимизирована работа с мультирастрами (не часто пока требуются).
В остальном инструмент меня устраивает, кого заинтересовал — отвечу на возникшие вопросы.