«Не все панки пришли к единому мнению о том, как нужно поддерживать других или каким образом можно что-то изменить во внешнем мире к лучшему. Но все согласны, что это необходимо делать.»
― Крейг О'Хара, «Философия панка: больше, чем шум»
Тему для второй заметки в рубрику «очумелые ручки» долго выбирать не пришлось — после прочтения первой, от одного из бета-тестеров нашего блога поступил вопрос: «А как в Python можно получить набор изолиний из грида для последующего использования их (изолиний) в геоинформационных системах по типу ArcGIS QGIS? Очень нада!!!»
Наш бродяга, вроде как загорелся идеей автоматизации процесса оформления отчетной графики в QGIS. Что ж, похвально! Как говорится, чем меньше времени мы тратим на рутину, тем больше времени остается на секс, пиво и рок-н-ролл творческую работу и эксперименты. Тем более, что тема действительно актуальная, в связи с последними требованиями к предоставляемым отчетным материалам в РГФ. Так давайте же убьем сразу двух зайцев: поможем братику и попутно разберем некоторые дополнительные типы пространственных данных, о которых мы не успели упомянуть в прошлый раз.
Let’s get started! Подгрузим для начала нашу поверхность, чтобы было с чем работать. Это мы с вами уже умеем. В качестве подопытного образца будем использовать карту изохрон из предыдущего набора данных:
|
|
Визуализировать двумерный массив данных (грид) с помощью библиотеки Matplotlib можно несколькими способами:
- Визуализация двумерного массива данных в виде растрового изображения — функция imshow библиотеки
Matplotlibили встроенная в библиотекуXarrayфункция plot (под капотом она тоже использует функциюimshowбиблиотекиMatplotlib) - Визуализация карты в виде контурных графиков с заливкой по заданным интервалам (векторная графика) с помощью функции contourf
- Визуализация изолиний — функция contour

Комбинирование различных вариантов визуализации на одном графике только приветствуется. Совместим контурный график и изолинии:
|
|

На вход функциям contour и contourf, необходимо передать одномерный массив уровней levels. Именно он определяет шаг между изолиниями и уровни закраски для контурного графика. Я создал его с помощью функции arange из библиотеки numpy, задав равномерный шаг в 10 миллисекунд для интервала между граничными значениями массива от -1700 до -1900 мс. Граничные значения уровней я задал вручную, хотя лучше было бы получить их непосредственно из грида, например вот так:
|
|
Окай, изолинии и векторную заливку мы отрисовали. И тут может возникнуть закономерный вопрос: так а если Matplotlib уже умеет извлекать изолинии из 2D массива данных, то не можем ли мы извлечь отрисованное для собственных нужд? Можем. Именно для этого мы предусмотрительно сохранили «слои» с изолиниями и контурной заливкой в переменные isoline и collec_poly. А теперь давайте посмотрим, что там внутри?
LineСтрингиString. Разбираем карту на изолинии
На первый взгляд, библиотека Matplotlib имеет довольно сложную структуру, а один вид ее документации может заставить вас рыдать в подушку, как девчонку. Читать документацию от корки до корки мы с вами конечно же не будем. Вместо этого попробуем применить грубую мужскую силу и решить нашу задачу методом реверс инжиниринга научного «тыка», попутно обращаясь к документации библиотеки (только в случае острой необходимости).
Вызовем переменную isoline и посмотрим, что нам вернет интерпретатор Python:
<matplotlib.contour.QuadContourSet at 0x1e40f8f2888>
Ничего не понятно, но очень интересно… Копируем данную строку в поисковик, чтобы получить больше информации. Получаем ответ: matplotlib.contour.QuadContourSet — это класс в Matplotlib, который представляет собой набор контурных линий или заполненных областей, сгенерированных функциями contour() или contourf(), а именно их мы использовали для создания контурных графиков и изолиний из 2D-данных.
Иными словами, как только мы вызываем функцию contour() или contourf(), она возвращает нам экземпляр класса QuadContourSet, который инкапсулирует (что?! О_о) всю информацию о сгенерированных контурах. Если еще проще: QuadContourSet — это грёбаная коробка, в которой хранится информация о контурных графиках, из которой в любой момент можно достать необходимое для дальнейших манипуляций. Именно то, что нам нужно.
Справедливости ради нужно отметить, что коробка эта имеет дополнительные отделения — атрибуты, в которых хранится разная информация о наших контурных графиках. С атрибутами класса QuadContourSet можно ознакомиться здесь. Нам пригодятся эти:
- сollections: список объектов
LineCollection(дляcontour()) илиPolyCollection(дляcontourf()), которые представляют собой фактические графические элементы контуров - levels: список числовых значений уровней контуров\изолиний
- allsegs: список, содержащий координаты сегментов контура для каждого уровня
Теперь, когда у нас есть больше информации, давайте попробуем вытащить что-нибудь из этого пакета с пакетами.
|
|
Данный фрагмент когда вернет нам следующее:
<a list of 20 PathCollection objects>
Python вернул нам список из двадцати PathCollection. Прям кроличья нора какая-то… Но нас это не остановит! Давайте двигаться последовательно: это список, а значит мы может получить конкретный его элемент по индексу. Давайте извлечем нулевой элемента списка:
|
|
Получаем:
<matplotlib.collections.PathCollection at 0x1e40f90a848>
Осталось выяснить что такое PathCollection. Идем в документацию, читаем. PathCollection — это набор (коллекция) линейных объектов Path. Да ладно! Они могут быть линейными, криволинейными, замкнутыми, разъединенными. Базовое хранилище для этих объектов включает в себя два numpy массива данных:
- vertices: в нем хранятся координаты вершин линейных объектов (в нашем случае изолиний)
- codes:
numpyмассив, который содержит в себе вспомогательную информацию, которая используетсяMatplotlibдля отрисовки линейных объектов, в зависимости от их типа. Что характерно, оба массива имеют одинаковую размерность.
Последнее с чем осталось разобраться, прежде чем мы наконец-то перейдем к коду: слово Collection в названии PathCollection как бы намекает нам, что внутри него не один Path. You know? Дело в том, что изолинии внутри Collections группируются по значению изолиний — levels. Проверим:
|
|
|
|
Итого: 20 PathCollection (ровно по количеству уровней наших изолиний), каждая из которых может содержать несколько изолиний (Path) одного уровня. Значит для того, чтобы собрать все это дерьмо в одну кучу нам потребуется запихнуть цикл в цикл. По итогу: один цикл будем проходиться по каждой из 20 PathCollection, а вложенный цикл будет собирать информацию о каждом Path. Фух!
Кстати, получить те самые два numpy массива vertices и codes из PathCollection можно с помощью функции get_paths(). Изолинии будем хранить в GeoDataFrame в виде геометрического примитива линии LineString. Теперь точно все. Погнали!
|
|
Ахтунг!
Должен предупредить: документация Matplotlib советует нам не обращаться к массивам вершин и кодов напрямую. Вместо этого, рекомендует обращаться к обозначенным массивам через операторы
iter_segmentsилиcleaned, которые помогают избежать проблем с None значениями в некоторых случаях. Но когда настоящих панков останавливали правила? Я, честно говоря, не вижу веских причин для того, чтобы придерживаться этого совета. Просто будем держать эту информацию в уме, а когда придет время воспользоваться этими знаниями, всенепременно так и поступим. Модификация кода, в случае необходимости, не потребует особых усилий.

Ну вот в общем-то и все. Можно сохранять в shp-файл. Хотя… Вы же не думаете, что мы остановимся на этом? Это было бы слишком просто для нас. На кой черт нам эти бездушные изолинии?! Давайте немного усложним задачу. Было бы неплохо, если бы наш скрипт дополнительно искал замкнутые контуры изолиний и сохранял их в отдельный полигональный GeoDataFrame, да еще и площадь этих структур считал.
С чего начнем? Для начала нужно понять, как мы будем определять замкнута ли изолиния или нет. Один из вариантов — использовать информацию из массива codes. Тем более, что сам Matplotlib использует его для отрисовки линий. Условно, по номеру кода в массиве codes Matplotlib понимает как соединять между собой две соседние точки из массива вершин vertices. Если пораскинуть мозгами, то не сложно догадаться, что нас интересует код последней точки массива codes у каждого Path. И тут все просто: если массив codes заканчивается цифрой ====79====, то линия считается замкнутой. Почему именно ====79====, спросите вы меня? Йух знает! Так написано в документации Matplotlib. Подробности можно найти тут.
По идее, все что нужно сделать — это добавить условие во вложенный цикл for с проверкой последнего значения в массиве codes для каждого линейного объекта Path. Если последнее значение равно ====79==== — считаем, что контур замкнут, в противном случае нет. Результат проверки пока будем сохранять в отдельный столбец нашего GeoDataFrame с изолиниями.
|
|
Посмотрим что получилось. Визуализируем наши данные, раскрасив изолинии в соответствии со значением из поля type:
|
|

Так, стоп. А это чей валет из колоды выпал?

И он там не один такой, если присмотреться внимательно. Первое, что приходит в голову: наверняка в этом контуре есть разрыв. В ходе внимательного визуального анализа контура — разрывов обнаружено не было. Кстати, для того, чтобы запустить график в интерактивном режиме (с возможностью увеличения необходимых областей и т.д.) необходимо использовать магическую команду %matplotlib notebook. Просто добавьте ее в ячейку с кодом Jupyter Notebook.
Вернемся к нашей проблеме: давайте посмотрим, что содержится в массивах vertices и codes замкнутого контура, тип которого определяется как «открытый». Я взял маленький контур из 5 точек. В нашем массиве GeoDataFrame он значится под индексом 39:
|
|

В структуре Matplotlib он находится в коллекции (PathCollection) с индексом 4 (изолинии со значением отметок -1850 м.), полигон с индексом 13 (Path):
|
|
|
|
Видно, что массив codes заканчивается кодовым значением 2, и по логике Matplotlib не является замкнутой линией. При этом координаты первой и последней точек (массив vertices) совпадают вплоть до последнего знака, что говорит о том, что данная линия является замкнутой. То есть она замкнута и не замкнута одновременно… Прям линия Шредингера какая-то.
Ок, пробежимся по всему пройденному «маршруту» координат данной линии в нашем коде. Посмотрим в промежуточном списке, куда мы сохраняли наши Path из цикла (переменная line_coll):
|
|
|
|
Все ок. Осталось посмотреть координаты этой записи в объекте GeoDataFrame:
|
|
|
|
Какого черта?! В каждой второй координате X и Y появилось три лишних знака после запятой (было 8, а стало 11)! И ТЕПЕРЬ уже координаты X первой и последней точек совпадают, а координаты Y этих же точек ОТЛИЧАЮТСЯ, причем именно в последних двух знаках после запятой. Тех самых, которые появляются как по мановению волшебной палочки, ёпта. Да и вообще, мне одному кажется странным, что во всех массивах координат отличается количество знаков после запятой через одну вершину?!!
Грёбаный Path этого Matplotlib!!!! И как с этим жить? Как по мне, так тут два варианта: либо мы с вами поймали bug-а за яйца, либо где-то жидко обосрались испражнились под себя и нужно было все-таки читать документацию более внимательно. В любом случае, отступать мы не собираемся! Предлагаю просто переписать часть кода с условием: теперь вместо того, чтобы проверять последнее значение из массива codes мы будем сравнивать координаты первой и последней точек из массива vertices. Если координаты будут совпадать — будем считать, что линия является замкнутой. А для пущей надежности немного округлим все значения координат, скажем до третьего знака после запятой.
Polygon. Определяем тип замкнутого контура
Для удобства, замкнутые изолинии будем собирать в отдельный полигональный GeoDataFrame. Для описания геометрии замкнутых полигональных объектов в GeoPandas используется примитив Polygon. Объекты Polygon создаются по аналогии с LineString, так что ничего принципиально нового тут не будет. Пока.
|
|

Ну вот, совсем другое дело! Все валеты в колоде. И теперь у нас есть отдельный GeoDataFrame с замкнутыми изолиниями — closed_gdf. Осталось научиться отличать положительные замкнутые (антиклинальные) структуры от отрицательных. Самая простая идея, которая пришла мне в голову — вычислить центральную точку (центроид) для каждого полигона и сравнивать значение отметок в этой точке со значением самого контура. Давайте по шагам:
- Вычисляем центроид для каждого полигона
- Снимаем значение абсолютной отметки с карты в точке центроида (этому мы научились в прошлой заметке)
- Сравниваем значение изолинии со значением абсолютной отметки центроида: если значение абсолютной отметки центроида больше, чем значение изолинии замкнутого контура, то структура считается положительной (антиклинальной). В противном случае — отрицательной.
Для вычисления центроида по геометрии объекта в GeoPandas есть стандартная функция centroid. Но есть нюанс: использование centroid не гарантирует нам, что точка попадет во внутреннюю область нашего полигона. Представьте замкнутый контур в виде подковы: вычисление centroid для такого полигона приведет к тому, что его центральная точка окажется вне контура полигона, и попадет ему в аккурат промеж булок. Для того, чтобы центральная точка гарантировано оказалась внутри контура полигона, в GeoPandas есть альтернативная функция representative_point().
Для вычисления площади полигонов в GeoPandas также есть специальная функция — area. Вот так просто. В этом и заключается все преимущество пространственного типа данных. Теперь мы знаем все необходимое, чтобы закончить наш скрипт:
|
|

Клац-клац и вот у нас: краткая сводка по структурам, необходимые shp-файлы для оформления графики в ГИС-проекте и несколько часов сэкономленного времени.

Сохраним shp-файлы
|
|
MultiPolygon. Векторная заливка карты
А что насчет цветной заливки? Да все аналогично! Без лишних слов напишем функцию, которая на вход будет принимать контурный график Matplotlib, а на выходе выплёвывать полигональный shp-файл с векторной заливкой нашей карты:
|
|
В отдельное поле мы сохранили цвет каждой области в hex формате, использовав функцию get_facecolor(). А геометрию полигонов, в этот раз, хранили в сгруппированном виде. В GeoPandas для этого есть специальный тип геометрии MultiPolygon. Аналогично работает MultiLineString для линейных объектов. Итак, посмотрим, что получили:
|
|

Итак, сегодня мы разобрали, как превратить грид в изолинии и подготовить их для работы в QGIS — теперь наш бета сигма-тестер сможет тратить меньше времени на рутинную графику и больше на действительно важные вещи. А заодно мы лишний раз убедились: Python и геоданные — это мощный дуэт для автоматизации, который не только экономит нервы, но и открывает пространство для творчества. Так что дерзайте, экспериментируйте, а если упретесь в новые вопросы — пишите. Ведь, как гласит панк-максима из эпиграфа: «что-то менять к лучшему — необходимо», даже если это начинается с таких «очумелых» рутинных задач.