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

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Нет описания правки
Строка 5: Строка 5:
== Знакомство с raster_tools ==
== Знакомство с raster_tools ==


В течении продолжительного времени развиваю инструментарий для препарирования георастров. Который является набором рецептов упрощающих взаимодействие с гибкой но сложной структурой python gdal. К сегодняшнему времени rater_tools напоминает некий «швейцарский нож» с разнообразным функционалом.  
В течении продолжительного времени развиваю инструментарий для препарирования георастров. Который является набором рецептов, упрощающих взаимодействие с гибкой, но сложной структурой Python + GDAL. К сегодняшнему времени '''raster_tools''' напоминает некий «швейцарский нож» с разнообразным функционалом.  
Код проекта расположен на [https://github.com/oldbay/raster_tools https://github.com/oldbay/raster_tools] , самый простой способ его установки в каталоге пользователя - при помощи pip:
Код проекта расположен на [https://github.com/oldbay/raster_tools https://github.com/oldbay/raster_tools] , самый простой способ его установки в каталоге пользователя - при помощи pip:


Строка 14: Строка 14:
'''raster_tools''' состоит из следующих классов:
'''raster_tools''' состоит из следующих классов:


'''raster2array''' - Базовый класс, отвечающий за весь функционал загрузки растра и сохранение numpy массива выбранного канала(band). Выгруженный из растра массив может быть преобразован методами данного класса.  
'''raster2array''' - Базовый класс, отвечающий за весь функционал загрузки растра и сохранение numpy-массива выбранного канала(band). Выгруженный из растра массив может быть преобразован методами данного класса.  


'''array2raster(raster2array)''' - Субкласс, отвечающий за сохранение обработанного массива в файле георастра. Он так же может сформировать виртуальный растр в памяти с возможностью обработки методами базового класса.
'''array2raster(raster2array)''' - Субкласс, отвечающий за сохранение обработанного массива в файле георастра. Он также может сформировать виртуальный растр в памяти с возможностью обработки методами базового класса.


'''raster2transform(raster2array)''' - Субкласс, отвечающий за трансформацию растра: изменение размерности массива и перепроецирование. Так же есть поддержка методов базового класса.
'''raster2transform(raster2array)''' - Субкласс, отвечающий за трансформацию растра: изменение размерности массива и перепроецирование. Также есть поддержка методов базового класса.


'''raster2calc(raster2array)''' - Класс итерационного калькулятора. Он не наследует базовый класс raster2array, но использует некоторые его методы опосредованно.
'''raster2calc(raster2array)''' - Класс итерационного калькулятора. Он не наследует базовый класс raster2array, но использует некоторые его методы опосредованно.
Строка 26: Строка 26:
'''multiarray2multiraster''' - Класс сохранения массивов, обработанных методами raster2multiarray, в мультиканальный растр.
'''multiarray2multiraster''' - Класс сохранения массивов, обработанных методами raster2multiarray, в мультиканальный растр.


Основной смысл существования всего этого «зоопарка» классов - это загрузка, обработка и выгрузка двухмерного numpy массива георастров. Помимо объекта 2-х мерного numpy массива некоторые методы raster2array позволяют выгружать данные в формате самоназвоного «стандартного словаря», имеющего следующий вид:
Основной смысл существования всего этого «зоопарка» классов - это загрузка, обработка и выгрузка двухмерного numpy-массива георастров. Помимо объекта 2-х мерного numpy-массива некоторые методы raster2array позволяют выгружать данные в формате самоназванного «стандартного словаря», имеющего следующий вид:


<syntaxhighlight lang="python">
<syntaxhighlight lang="python">
Строка 39: Строка 39:
Эта конструкция несёт в себе всю необходимую информацию для создания нового георастра, отличного от исходного.
Эта конструкция несёт в себе всю необходимую информацию для создания нового георастра, отличного от исходного.


Углубляться в подробное описание стрктуры классов и методов raster_tools не вижу необходимости - достаточно полная базовая документация есть в README.rst на [https://github.com/oldbay/raster_tools https://github.com/oldbay/raster_tools] . В данной статье хотелось бы показать «боевое» применение raster_tools на практике, в том виде как сам его в проектах «готовлю».
Углубляться в подробное описание структуры классов и методов raster_tools нет необходимости - достаточно полная базовая документация есть в README.rst на [https://github.com/oldbay/raster_tools https://github.com/oldbay/raster_tools]. В данной статье хотелось бы показать «боевое» применение raster_tools на практике.


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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 47: Строка 47:
</syntaxhighlight>
</syntaxhighlight>


В ветке по умолчанию (master) , расположены только тестовые скрипты и исходные данные. Чтобы заранее изучить результаты тестов, нужно перейти в ветку result:
В ветке по умолчанию (master) расположены только тестовые скрипты и исходные данные. Чтобы заранее изучить результаты тестов, нужно перейти в ветку "result":


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 53: Строка 53:
</syntaxhighlight>
</syntaxhighlight>


Тестовый репозирорий состоит из:
Тестовый репозиторий состоит из:
* Корень: содержит тестовые примеры - выполнение которых необходимо производить из текущего каталога;
* Корень: содержит тестовые примеры, выполнение которых необходимо производить из текущего каталога;
* Каталог data: содержит исходные данные для экспериментов;
* Каталог data: содержит исходные данные для экспериментов;
* Каталог logs: содержит результаты тестов расхода памяти и времени выполнения (branch result);
* Каталог logs: содержит результаты тестов расхода памяти и времени выполнения (branch result);
* Каталог module: содержит дополнительные python модули, необходимые для выполнение некоторых тестовых скриптов(в статье не описаны);
* Каталог module: содержит дополнительные Python-модули, необходимые для выполнение некоторых тестовых скриптов(в статье не описаны);
* Каталог results: содержит результаты экспериментов (branch result).
* Каталог results: содержит результаты экспериментов (branch result).


Строка 65: Строка 65:
Первоначально raster_tools был частью растрового калькулятора, сейчас отдельный проект raster_calc, поэтому исторически начнём с вычислений.
Первоначально raster_tools был частью растрового калькулятора, сейчас отдельный проект raster_calc, поэтому исторически начнём с вычислений.
   
   
Пример простого растрового калькулятора представлен в calc.py (list 1) . В данном скрипте мультирастр разбирается на спектральные каналы, после чего вычисляется вегетативный индекс на основе rgb - TGI([https://gist.github.com/merkato/0cd894f19518496171afd7425e09ed88 украдено отсюда]), Источником данных является data/multi.tif (img 1).
Пример простого растрового калькулятора представлен в calc.py (list 1) . В данном скрипте мультирастр разбирается на спектральные каналы, после чего вычисляется вегетативный индекс на основе rgb - TGI ([https://gist.github.com/merkato/0cd894f19518496171afd7425e09ed88 украдено отсюда]), Источником данных является data/multi.tif (img 1).




Строка 119: Строка 119:




У метода есть существенный минус - высокое потребление памяти. Фактически в памяти , на пике потребления, присутствуют сразу 4 массива растров: каналов r, g, b и вычисляемого индекса. При работе с большими георастрами это часто неприемлемо. Поэтому был разработан итерационный метод работы растрового калькулятора. Данный метод работает с растром посекторно в цикле. При таких вычислениях сильно снижается потребление памяти, но падает производительность.
У метода есть существенный минус - высокое потребление памяти. Фактически в памяти на пике потребления присутствуют сразу 4 массива растров: каналов r, g, b и вычисляемого индекса. При работе с большими георастрами это часто неприемлемо. Поэтому был разработан итерационный метод работы растрового калькулятора. Данный метод работает с растром посекторно в цикле. При таких вычислениях сильно снижается потребление памяти, но падает производительность.


Итерационный калькулятор реализован в calc_iter.py (list 2). В отличии от простого калькулятора здесь формула вычислений(в формате lambda) и переменные для вычислений(в виде объектов raster2array) передаются классу raster2calc.
Итерационный калькулятор реализован в calc_iter.py (list 2). В отличии от простого калькулятора здесь формула вычислений(в формате lambda) и переменные для вычислений(в виде объектов raster2array) передаются классу raster2calc.
Строка 163: Строка 163:
     )
     )


# сохраниение стандартного словаря в растр.
# сохранение стандартного словаря в растр.
array2raster(None, out, out_file)
array2raster(None, out, out_file)
</syntaxhighlight>
</syntaxhighlight>


raster2calc возвращает «стандартный словарь» преобразуемый классом array2raster в файл result/calc_iter.tif.
raster2calc возвращает «стандартный словарь», преобразуемый классом array2raster в файл result/calc_iter.tif.


Если выполнить calc.py и calc_iter.py в обёртке memory_profiler - получаем графики потребления памяти:
Если выполнить calc.py и calc_iter.py в обёртке memory_profiler - получаем графики потребления памяти:
Строка 174: Строка 174:
[[Файл:Calc_iter_mprof.png|600px]]
[[Файл:Calc_iter_mprof.png|600px]]


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


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


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


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


Стандартный словарь применим и в растровом калькуляторе, для подобных вычислений необходимо вернуть numpy массив по ключу «array» list 1.1
Стандартный словарь применим и в растровом калькуляторе, для подобных вычислений необходимо вернуть numpy массив по ключу «array» list 1.1
Строка 197: Строка 197:
blue = raster2array(in_file, 3).get_std_dict()
blue = raster2array(in_file, 3).get_std_dict()


# сохранение numpy массивов каналов в переменные.
# сохранение numpy-массивов каналов в переменные.
r = red()["array"]
r = red()["array"]
g = green()["array"]
g = green()["array"]
Строка 207: Строка 207:




# сохранение вычисленного массива в стандартный словарь
# сохранение вычисленного массива в стандартный словарь:
out_std_dict = red
out_std_dict = red
out_std_dict["array"] = calc
out_std_dict["array"] = calc


# сохранение стандартного словаря в растр.
# сохранение стандартного словаря в растр:
array2raster(None, out_std_dict, out_file)
array2raster(None, out_std_dict, out_file)
</syntaxhighlight>
</syntaxhighlight>


Теперь рассмотрим работу каждого метода обрезки растра:
Теперь рассмотрим работу каждого метода обрезки растра.




=== Метод cut_area. ===
=== Метод cut_area ===


Позволяет обрезать растр по набору координат точек. Применяется когда необходимо обрезать область растра по массиву точек - для дальнейшей обработки. Пример использования приведён в cut_coords.py (list 3), в качестве исходных данных используется result/calc.tif.
Позволяет обрезать растр по набору координат точек. Применяется когда необходимо обрезать область растра по массиву точек - для дальнейшей обработки. Пример использования приведён в cut_coords.py (list 3), в качестве исходных данных используется result/calc.tif.
Строка 270: Строка 270:
out_file = "result/cut_shp.tif"
out_file = "result/cut_shp.tif"


# вырезание области растра по полигону shp файла
# вырезание области растра по полигону shp-файла
# и сохранение в формате стандартного словаря
# и сохранение в формате стандартного словаря
std_dict = raster2array(in_file).cut_shp_file(shp_file)
std_dict = raster2array(in_file).cut_shp_file(shp_file)


# сохранеие стандартного словаря в растр
# сохранение стандартного словаря в растр
array2raster(None, std_dict, out_file)
array2raster(None, std_dict, out_file)
</syntaxhighlight>
</syntaxhighlight>
Строка 281: Строка 281:




img 4.1 - геометрия из shp файла на фоне исходного растра
img 4.1 - геометрия из shp-файла на фоне исходного растра


[[Файл:Cut_shp_screenshot.png|300px]]
[[Файл:Cut_shp_screenshot.png|300px]]
Строка 291: Строка 291:




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




=== Метод cut_ogr_geometry. ===
=== Метод cut_ogr_geometry ===


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


Для выполнения примера следует провести подготовку(все примеры работы с БД приведены для linux):  
Для выполнения примера следует провести подготовку(все примеры работы с БД приведены для Linux):  
* Установить postgresql и расширение postgis.  
* Установить PostgreSQL c расширением PostGIS.  
* В СУБД создать роль из под которой возможен логин в БД для загрузки дампа:
* В СУБД создать роль, из-под которой возможен логин в БД для загрузки дампа:
<pre>
<pre>
# su - postgres
# su - postgres
Строка 319: Строка 319:
psql -d doc_example < data/crowns.sql
psql -d doc_example < data/crowns.sql
</pre>
</pre>
* Проверить наличие python модуля psycopg2 - драйвера postgresql.(необходим для интерфейса запросов modules.db_interface)
* Проверить наличие Python-модуля "psycopg2" - драйвера PostgreSQL (необходим для интерфейса запросов modules.db_interface)


После успешного завершения подготовительного этапа можно выполнять тестовый пример cut_wkt.py (list 5)
После успешного завершения подготовительного этапа можно выполнять тестовый пример cut_wkt.py (list 5)
Строка 345: Строка 345:
geom_id = 55775
geom_id = 55775


# запрос в БД возвращающий полигон в формате WKT
# запрос в БД, возвращающий полигон в формате WKT
SQL = """
SQL = """
select ST_AsText(wkb_geometry)
select ST_AsText(wkb_geometry)
Строка 385: Строка 385:




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


Иногда приходится вычислять индексы, или иным способом сравнивать георастры из различных источников. Входные данные могут быть в разном масштабе и соответственно иметь разную размерность массивов исследуемой области. В результате вычисление индексов средствами numpy становится невозможным, а многие другие операции затруднительны.  
Иногда приходится вычислять индексы, или иным способом сравнивать георастры из различных источников. Входные данные могут быть в разном масштабе и соответственно иметь разную размерность массивов исследуемой области. В результате вычисление индексов средствами numpy становится невозможным, а многие другие операции затруднительны.  
Строка 444: Строка 444:
[[Файл:Resize_300_300.png|500px]]
[[Файл:Resize_300_300.png|500px]]


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


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


* raster2multiarray - Класс умеющий изменять очерёдность и количество каналов, при этом преобразующий многоканальный растр в стандартный словарь расширенного типа. От обычного он отличается тем что в поле «array» сохраняется не двухмерный (x, y), а трёхмерный numpy массив с индексом (bands, x, y). Для определения типа массива добавлено поле «multi_type» - принимающее значения: None - с массивом уже указанного индекса и «cv» - с индексом (x, y, bands). Массив типа multi_type = «cv» аналогичен возвращаемому функцией cv2.imread. Кроме того у raster2multiarray имеется набор методов обрезки растра аналогичных raster2array.
* raster2multiarray - Класс, умеющий изменять очерёдность и количество каналов, при этом преобразующий многоканальный растр в стандартный словарь расширенного типа. От обычного он отличается тем, что в поле «array» сохраняется не двухмерный (x, y), а трёхмерный numpy-массив с индексом (bands, x, y). Для определения типа массива добавлено поле «multi_type», принимающее значения: None - с массивом уже указанного индекса и «cv» - с индексом (x, y, bands). Массив типа multi_type = «cv» аналогичен возвращаемому функцией cv2.imread. Кроме того, у raster2multiarray имеется набор методов обрезки растра, аналогичных raster2array.


* multiarray2multiraster - Класс сохраняющий многоканальный стандартный словарь в мультиканальный растр.
* multiarray2multiraster - Класс, сохраняющий многоканальный стандартный словарь в мультиканальный растр.


Предлагаю два варианта использования данной пары классов, исходным растром для примеров послужит data/multi.tif:
Предлагаю два варианта использования данной пары классов, исходным растром для примеров послужит data/multi.tif:




===Изменение последовательности каналов.===
=== Изменение последовательности каналов ===


Для данного примера в rgb мультирастре переставим местами 1-й и 3-й каналы. Тестовый скрипт представлен в multi_rebands.py (list 7).  
Для данного примера в RGB-мультирастре переставим местами 1-й и 3-й каналы. Тестовый скрипт представлен в multi_rebands.py (list 7).  


list 7 - multi_rebands.py
list 7 - multi_rebands.py
Строка 473: Строка 473:
img = raster2multiarray(img_in, 3, 2, 1)
img = raster2multiarray(img_in, 3, 2, 1)


# выгрузка многоканального стандартого словаря
# выгрузка многоканального стандартного словаря
img = img.get_std_dict()
img = img.get_std_dict()


# cохранение многоканального стандартого словаря
# cохранение многоканального стандартного словаря
# в мультирастр
# в мультирастр
multiarray2multiraster(img_out, img)
multiarray2multiraster(img_out, img)
Строка 489: Строка 489:




===Обработка муьтирастров в opencv.===
=== Обработка муьтирастров в opencv ===


Теперь усложним задачу: удалим альфа канал мультирастра, передадим массив на обработку opencv и сохраним полученный результат. Для этого примера необходимо установить python opencv 2 или 3. Тестовый скрипт представлен в multi_opencv.py (list 8).
Теперь усложним задачу: удалим альфа-канал мультирастра, передадим массив на обработку OpenCV и сохраним полученный результат. Для этого примера необходимо установить Python-модули OpenCV версии 2 или 3. Тестовый скрипт представлен в multi_opencv.py (list 8).


list 8 - multi_opencv.py
list 8 - multi_opencv.py
Строка 539: Строка 539:
==Работа с точками.==
==Работа с точками.==


Не очевидный, но временами очень важный, метод препарирования георастров. Постановка задачи: запись значений точек растра с определёнными координатами в csv, БД или в виде поля в слой shp файла. Реализация последнего примера представлена в point_data_load_array.py (list 9). Исходными данными здесь выступают: георастр  result/calc.tif и векторный файл со слоем точек data/tops.shp.
Неочевидный, но временами очень важный метод препарирования георастров. Постановка задачи: запись значений точек растра с определёнными координатами в CSV, БД или в виде поля в слой shp-файла. Реализация последнего примера представлена в point_data_load_array.py (list 9). Исходными данными здесь выступают: георастр  result/calc.tif и векторный файл со слоем точек data/tops.shp.


list 9 - point_data_load_array.py
list 9 - point_data_load_array.py
Строка 552: Строка 552:
raster = "result/calc.tif"
raster = "result/calc.tif"


# копируем тестовый shp файл с точками в каталог result
# копируем тестовый shp-файл с точками в каталог result
_filename = "tops"
_filename = "tops"
_indir = "data"
_indir = "data"
Строка 570: Строка 570:
raster.np_array_load()
raster.np_array_load()


# загрузка shp файла
# загрузка shp-файла
_shp = ogr.Open(shp_file, update=1)
_shp = ogr.Open(shp_file, update=1)
layer = _shp.GetLayerByIndex(0)
layer = _shp.GetLayerByIndex(0)


# добавление в shp файл нового поля
# добавление в shp-файл нового поля
# для записи выгруженных значений растра
# для записи выгруженных значений растра
new_field = "POINT_DATA"
new_field = "POINT_DATA"
Строка 581: Строка 581:
     layer.CreateField(field)
     layer.CreateField(field)


# обход циклом точек shp файла
# обход циклом точек shp-файла
for geometry in layer:
for geometry in layer:


Строка 592: Строка 592:


     # сохранение выгруженных из растра данных
     # сохранение выгруженных из растра данных
     # в новое поле shp файла
     # в новое поле shp-файла
     geometry.SetField(new_field, raster_data)
     geometry.SetField(new_field, raster_data)
     layer.SetFeature(geometry)
     layer.SetFeature(geometry)
Строка 606: Строка 606:
</syntaxhighlight>
</syntaxhighlight>


Такая конструкция загружает numpy массив растра непосредственно в экземпляр класса rater2array - что уменьшает время выполнения операций, увеличивая потребление памяти. Для сравнения: в репозитории есть аналогичный скрипт point_data_dont_load_array.py где данная строка отсутствует. Сравнение графиков выполнения обоих скриптов в memory_profiler выглядит следующим образом:  
Такая конструкция загружает numpy-массив растра непосредственно в экземпляр класса raster2array, что уменьшает время выполнения операций, увеличивая потребление памяти. Для сравнения: в репозитории есть аналогичный скрипт point_data_dont_load_array.py, где данная строка отсутствует. Сравнение графиков выполнения обоих скриптов в memory_profiler выглядит следующим образом:  


[[Файл:Point_data_load_array_mprof.png|600px]]
[[Файл:Point_data_load_array_mprof.png|600px]]
Строка 612: Строка 612:




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


Иногда процедурная генерация георастров некоторыми утилитами порождает «мутантов». При их отображении в gis программах проблема не заметна - так как геоданные будут правильно спозицианированны. Но при выгрузке из такого георастра numpy массива - оказывается что он зеркально отображён относительно оси x или y. Примером утилиты с такой странной процедурной генерации является gdal_grid([https://gis-lab.info/forum/viewtopic.php?f=30&t=20283 сообщение о проблеме на форуме]).
Иногда процедурная генерация георастров некоторыми утилитами порождает «мутантов». При их отображении в ПО ГИС проблема не заметна - так как геоданные будут правильно позиционированы. Но при выгрузке из такого георастра numpy-массива оказывается, что он зеркально отображён относительно оси X или Y. Примером утилиты с такой странной процедурной генерации является gdal_grid ([https://gis-lab.info/forum/viewtopic.php?f=30&t=20283 сообщение о проблеме на форуме]).
   
   
Чтобы не быть голословными, выполним весь цикл процедур: интерполяция данных из предыдущего примера result/tops.shp при помощи gdal_grid, проверка полученного процедурного растра на «зеркальность», исправление его при помощи raster_tools. В качестве исходных данных будут использованы: result/tops.shp (источник данных для интерполяции) и result/calc.tif (для сравнения). Обработку произведём скриптом repair.py (list 11). Для работы repair.py необходим модуль modules.grid_utils и сама утилита gdal_rid.   
Чтобы не быть голословными, выполним весь цикл процедур: интерполяция данных из предыдущего примера 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
list 11 - repair.py
Строка 633: Строка 633:
repair_file = "result/repair.tif"
repair_file = "result/repair.tif"


# Интерполяция на основе shp файла
# Интерполяция на основе shp-файла
grid_utils(raster_file, shp_file, field_name, out_dir)
grid_utils(raster_file, shp_file, field_name, out_dir)


Строка 649: Строка 649:
</syntaxhighlight>
</syntaxhighlight>


В результате работы скрипта были созданы 2 растра: result/POINT_DATA.tif - результат работы gdal_grid и исправленный растр result/repair.tif. Сравним их gdalinfo лог:
В результате работы скрипта были созданы 2 растра: result/POINT_DATA.tif - результат работы gdal_grid и исправленный растр result/repair.tif. Сравним вывод информации о растрах через gdalinfo:


<pre>
<pre>
Строка 706: Строка 706:
</pre>
</pre>


Обратите внимание на второе значение Pixel Size(ось y): у растра сгенерированного gdal_grid это значение положительное, в отличии от исходного и исправленного георастров. Что и позволяет отображать POINT_DATA.tif в гис программах нормально, несмотря на наличие зеркально отображённого через ось y массива. То же самое сообщает нам вывод метода raster2array.is_valid():
Обратите внимание на второе значение Pixel Size(ось y): у растра, сгенерированного gdal_grid, это значение положительное в отличие от исходного и исправленного георастров. Что и позволяет отображать POINT_DATA.tif в ПО ГИС нормально, несмотря на наличие зеркально отображённого через ось Y массива. То же самое сообщает нам вывод метода raster2array.is_valid():


<syntaxhighlight lang="python">
<syntaxhighlight lang="python">
Строка 721: Строка 721:
</syntaxhighlight>
</syntaxhighlight>


Здесь ориентация растра по оси x истинное и ложное по оси y.
Здесь ориентация растра по оси X истинное и ложное по оси Y.
 


== Послесловие ==


==Послесловие.==
Заканчивая описание способов употребления raster_tools, хочется отметить, что инструмент будет и дальше находиться в состоянии постоянного развития. По мере поступления задач в него будут добавляться какие-то новые функции, при этом с сохранением старых методов применения (чтобы не потерять совместимость с уже завершёнными проектами). Статья описывает состояние raster_tools версии 0.4 - на текущий момент последняя стабильная версия.


Заканчивая описание способов употребления raster_tools хочу сказать что инструмент будет и дальше находится в состоянии постоянного развития. По мере поступления задач буду добавлять в него какие то новые приспособления, при этом постараюсь сохранить старые методы применения (чтобы не потерять совместимость с уже завершёнными проектами). Статья описывает состояние raster_tools версии 0.4 - на текущий момент последний stable.
У raster_tools имеются следующие недостатки:  
У raster_tools имеются следующие недостатки:  
* проблемное перепроецирование векторных слоёв-шаблонов в проекцию растра, поэтому стараюсь использовать единую проекцию;
* проблемное перепроецирование векторных слоёв-шаблонов в проекцию растра, поэтому логично использовать единую проекцию;
* обрезка по shp файлам «прибита гвоздями» к первому векторному слою;
* обрезка по shp-файлам «прибита гвоздями» к первому векторному слою;
* инструмент рассчитан на geotif - работа с другими растровыми форматами файлов в зачаточном состоянии;
* инструмент рассчитан на формат GeoTIFF; работа с другими растровыми форматами файлов в зачаточном состоянии;
* недостаточно оптимизирована работа с мультирастрами (не часто пока требуются).
* недостаточно оптимизирована работа с мультирастрами (нечасто пока требуются);
В остальном инструмент меня устраивает, кого заинтересовал — отвечу на возникшие вопросы.

Версия от 21:24, 30 июля 2018

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


В статье рассматриваются примеры практического использования набора растровых инструментов 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

Multi.png


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

Calc.png


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

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

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

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

Cut coords screenshot.png

img 3.2 - result/cut_coords.tif

Cut coords.png


Метод 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-файла на фоне исходного растра

Cut shp screenshot.png


img 4.2 - result/cut_shp.tif

Cut shp.png


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


Метод 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 на фоне исходного растра

Cut wkt screenshot.png


img 5.2 - result/cut_wkt.tif

Cut wkt.png


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

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

Resize 253 210.png


img 6.1 - разрешение (125x105)

Resize 125 105.png


img 6.1 - разрешение (300x300)

Resize 300 300.png

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

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

  • 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

Multi rebands.png


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

Multi opencv.png


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

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

Point data load array mprof.png Point data dont load array mprof.png


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

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