Практика использования raster tools: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
(Новая страница: «{{Статья|Черновик}} {{Аннотация|В статье рассматриваются примеры практического использов…»)
 
Нет описания правки
Строка 5: Строка 5:


В течении продолжительного времени развиваю инструментарий для препарирования георастров. Который является набором рецептов упрощающих взаимодействие с гибкой но сложной структурой python gdal. К сегодняшнему времени rater_tools напоминает некий «швейцарский нож» с разнообразным функционалом.  
В течении продолжительного времени развиваю инструментарий для препарирования георастров. Который является набором рецептов упрощающих взаимодействие с гибкой но сложной структурой python gdal. К сегодняшнему времени rater_tools напоминает некий «швейцарский нож» с разнообразным функционалом.  
Код проекта расположен на [https://github.com/oldbay/raster_tools 'github'] , самый простой способ его установки в каталоге пользователя - при помощи pip:
Код проекта расположен на [https://github.com/oldbay/raster_tools https://github.com/oldbay/raster_tools] , самый простой способ его установки в каталоге пользователя - при помощи pip:
 
<pre>
pip2.7 install --user git+https://github.com/oldbay/raster_tools
</pre>
 
'''raster_tools''' состоит из следующих классов:
 
'''raster2array''' - Базовый класс, отвечающий за весь функционал загрузки растра и сохранение numpy массива выбранного канала(band). Выгруженный из растра массив может быть преобразован методами данного класса.
 
'''array2raster'''(raster2array) - Субкласс, отвечающий за сохранение обработанного массива в файле георастра. Он так же может сформировать виртуальный растр в памяти с возможностью обработки методами базового класса.
 
'''raster2transform'''(raster2array) - Субкласс, отвечающий за трансформацию растра: изменение размерности массива и перепроецирование. Так же есть поддержка методов базового класса.
 
'''raster2calc'''(raster2array) - Класс итерационного калькулятора. Он не наследует базовый класс raster2array, но использует некоторые его методы опосредованно.
 
'''raster2multiarray'''(raster2array) - Класс обработки мультирастров. Позволяет миксовать каналы мультиканального растра и обрабатывать его некоторыми методами raster2array.
 
'''multiarray2multiraster''' - Класс сохранения массивов, обработанных методами raster2multiarray, в мультиканальный растр.
 
Основной смысл существования всего этого «зоопарка» классов - это загрузка, обработка и выгрузка двухмерного numpy массива георастров. Помимо объекта 2-х мерного numpy массива некоторые методы raster2array позволяют выгружать данные в формате самоназвоного «стандартного словаря», имеющего следующий вид:
 
<syntaxhighlight lang="python">
{
    "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()
}
</syntaxhighlight> 
 
Эта конструкция несёт в себе всю необходимую информацию для создания нового георастра, отличного от исходного.
 
Углубляться в подробное описание стрктуры классов и методов raster_tools не вижу необходимости - достаточно полная базовая документация есть в README.rst на [https://github.com/oldbay/raster_tools https://github.com/oldbay/raster_tools] . В данной статье хотелось бы показать «боевое» применение raster_tools на практике, в том виде как сам его в проектах «готовлю».
 
Для этих экспериментов был создан gis репозиторий с тестовыми данными , загрузить которые можно командой:
 
<pre>
git clone  https://gitlab.com/oldbay/raster_tools_examples.git
</pre>
 
В ветке по умолчанию (master) , расположены только тестовые скрипты и исходные данные. Чтобы заранее изучить результаты тестов, нужно перейти в ветку result:
 
<pre>
git checkout origin/result
</pre>
 
Тестовый репозирорий состоит из:
* Корень: содержит тестовые примеры - выполнение которых необходимо производить из текущего каталога;
* Каталог data: содержит исходные данные для экспериментов;
* Каталог logs: содержит результаты тестов расхода памяти и времени выполнения (branch result);
* Каталог module: содержит дополнительные python модули, необходимые для выполнение некоторых тестовых скриптов(в статье не описаны);
* Каталог results: содержит результаты экспериментов (branch result).
 
== 2. Растровый калькулятор. ==
 
Первоначально raster_tools был частью растрового калькулятора, сейчас отдельный проект raster_calc, поэтому исторически начнём с вычислений.
Пример простого растрового калькулятора представлен в calc.py (list 1) . В данном скрипте мультирастр разбирается на спектральные каналы, после чего вычисляется вегетативный индекс на основе rgb - TGI([https://gist.github.com/merkato/0cd894f19518496171afd7425e09ed88 украдено отсюда]), Источником данных является data/multi.tif.
 
list 1 - calc.py
<syntaxhighlight lang="python">
#!/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)
</syntaxhighlight>
 
В результате выполнения calc.py создается растр result/calc.tif.
 
Img:
 
У метода есть существенный минус - высокое потребление памяти. Фактически в памяти , на пике потребления, присутствуют сразу 4 массива растров: каналов r, g, b и вычисляемого индекса. При работе с большими георастрами это часто неприемлемо. Поэтому был разработан итерационный метод работы растрового калькулятора. Данный метод работает с растром посекторно в цикле.  При таких вычислениях сильно снижается потребление памяти, но падает производительность.
 
Итерационный калькулятор реализован в calc_iter.py (list 2). В отличии от простого калькулятора здесь формула вычислений(в формате lambda) и переменные для вычислений(в виде объектов raster2array) передаются классу raster2calc.
 
list 2 - calc_iter.py
<syntaxhighlight lang="python">
#!/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)
</syntaxhighlight>
 
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
 
Видно что итерационный метод в несколько раз выигрывает в потреблении памяти, но теряет в скорости расчётов. Потери скорости не очень существенны на текущем примере - но при массовой обработке небольших георастров, бывает выгоднее пожертвовать памятью.

Версия от 15:21, 31 мая 2018

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


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

1. Знакомство с 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).

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

Первоначально 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

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