Практика использования 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.

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:

У метода есть существенный минус - высокое потребление памяти. Фактически в памяти , на пике потребления, присутствуют сразу 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 - получаем графики потребления памяти:

calc.py

img: logs/calc_mprof.png

calc_iter.py

img:logs/calc_iter_mprof.png

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

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

Иногда растровые вычисления необходимо производить только с некоторой выделенной областью георастра. Для «фигурного» вырезания областей, ограниченных координатами точек или 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:

Метод 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:

Обратите внимание что растр обрезан непосредственно по полигону, а не только по координатам 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:

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

Иногда приходится вычислять индексы, или иным способом сравнивать георастры из различных источников. Входные данные могут быть в разном масштабе и соответственно иметь разную размерность массивов исследуемой области. В результате вычисление индексов средствами 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), (125,105), (300,300). Последний массив сохранится в растр result/resize.tif.

Img:

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

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

• 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:

Обработка муьтирастров в 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:

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

Не очевидный, но временами очень важный, метод препарирования георастров. Постановка задачи: запись значений точек растра с определёнными координатами в 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 выглядит следующим образом:

point_data_load_array.py

img: logs/point_data_load_array_mprof.png

point_data_dont_load_array.py

img: logs/point_data_dont_load_array_mprof.png