/img/star (2).png.webp)
/img/news (2).png.webp)
/img/event.png.webp)
Время прочтения: 13 мин.
Итак, суть публикации заключена в знакомстве с темой отображения Земли на плоскости, видами картографических проекций и их нюансами, инструментом Cartopy, позволяющим представить какие-либо объекты и области на карте.
Изложу Вам план работы (действуем согласно заголовку поста):
- «Поговорим» об отображении поверхности Земли на плоскости, основах представления карт и данных, привязанных к координатам/областям на них; рассмотрим примеры визуализаций – отобразим карты и данные с помощью модуля.
- Проверим возможности модуля на 4-х тривиальных кейсах, с которыми я столкнулся на практике.
- В заключение кратко изложим, что из себя представляет объект поста (Cartopy), пути его применения в сфере бизнеса и финансов, имеет ли смысл использовать этот инструмент для визуализации данных.
Cartopy — это Python-пакет, который в связке с Matplotlib позволяет визуализировать пространственные данные на картах — легко, наглядно и, главное, точно. Ключевые особенности модуля: объектно-ориентированное определение проекций и способность преобразовывать точки, линии, векторы, многоугольники и изображения между ними (проекциями). Cartopy используется, например, в метеорологии, изучении поверхности Земли и околоземного космического пространства.
Приступим к установке Cartopy (предварительно рекомендую создать виртуальное окружение).
Это можно сделать с помощью conda:
conda install -c conda-forge cartopyЛибо через pip, установив перед этим зависимости, без которых библиотека не сможет работать: Shapely, pyshp, pyproj и GEOS. Подчеркну, что GEOS (libgeos) и geos (pip geos) для нас не одно и то же. После успешной инсталляции GEOS, остальные зависимости присоединяем через pip:
pip install shapely pyshp pyprojНу и сам Cartopy:
pip install cartopyДалее подключим используемые для работы с данными и визуализацией библиотеки:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import cartopy.feature as sfeatureТак как Cartopy построен поверх интерфейса Matplotlib, его синтаксис должен быть знаком и прост для понимания, структура диаграмм аналогичная.
Напомню, что в Matplotlib.pyplot диаграмма имеет иерархическую структуру. Основным элементом является Figure (рисунок – объект верхнего уровня), на нём расположены Axes (области рисования, где изображается пространство данных), Axis (оси координат, которые упорядочивают отображение данных) и Artists (элементы рисунка – заголовок, легенда, линия, текст, деления ticks, подписи к ним tick labels и т.д.).
Основные структурные элементы представлены на рисунке ниже. Более подробно изучить структуру Matplotlib можно с помощью документации.
Ранее мы импортировали crs и feature — два основных компонента пакета. Первый содержит проекции с системами отсчёта координат и функцию преобразования данных между ними. Map projection или картографическая проекция — математическое преобразование широт/долгот местоположений объектов с поверхности сферы или эллипсоида на плоскость. По типу проецирования Map projections можно разделить на азимутальные, конические, цилиндрические и иные (псевдо-проекции).
Основные типы проекций – азимутальная, коническая и цилиндрическая
Из-за разной топологии любой проекции присущи искажения — в площади, форме, углах, расстоянии или масштабе. На данный момент не существует универсальной проекции, которая бы сохраняла параметры поверхности Земли без искажений. Одни проекции сохраняют конкретное свойство, другие — балансируют между неточностями. Поэтому важно выбирать «правильную» проекцию, исходя из задачи. Она, как и любая диаграмма, при некорректной репрезентации может привести к тому, что данные, полученные в результате качественного анализа, будут интерпретированы неверно.
Пробежимся по основам модуля: как отобразить проекцию, какие проекции представлены в Cartopy. В дальнейшем все примеры и задачи будут связаны с визуализацией проекции, данных на ней. Проекцию Земли можно построить в пару строк. Выведем несколько экземпляров класса crs.
[Код]
```
# задание области рисования c размерами 9x4.5 дюймов
fig = plt.figure(figsize=(9, 4.5))
# выделение заголовка
fig.suptitle(‘Projections’, fontsize=20, fontweight=’bold’, color=’gray’, x=0.5, y=0.97)
projections = {‘Robinson’: ccrs.Robinson(), ‘AlbersEqualArea’: ccrs.AlbersEqualArea(), ‘AzimuthalEquidistant’: ccrs.AzimuthalEquidistant(),
‘Mercator’: ccrs.Mercator(), ‘InterruptedGoodeHomolosine’: ccrs.InterruptedGoodeHomolosine(), ‘LambertCylindrical’: ccrs.LambertCylindrical()}
# проход по словарю для отображения проекций вместе с наименованиями
for index, proj in enumerate(projections.items()):
ax = fig.add_subplot(2, 3, index + 1, projection=proj[1])
# отображение границ суши
ax.coastlines()
ax.set_title(proj[0])
plt.tight_layout()
plt.show()
```
Второй импортированный компонент — features — рисует точки, линии, полигоны, тем самым добавляет стиль, детали на картах.
Обдуманно в дальнейшем будем использовать две проекции: Plate-Carree и Orthographic, чтобы удовлетворить ценителей плоской Земли и «шарящих» т.к. это одни из основных используемых картографических отображений поверхности Земли. Plate-Carree сохраняет неизменным масштаб линий направлений, удобна для работы с географической системой координат, поэтому применятся в геоинформационных системах (де-факто стандарт в ГИС). Ортографическая проекция сохраняет центр карты без искажений, используется для отображения полярных регионов, вида Земли из космоса.
Для первой из них применим «фичи»: отметим сушу, воду, береговую линию, выделим пунктиром границы государств через параметр BORDERS, с помощью STATE отобразим территориальные единицы США (для других стран схожего параметра пока нет).
[Код]
```
plt.figure()
ax = plt.axes(projection=ccrs.Plate-Carree)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=’:’)
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.STATES)
```
Также просто изобразить рельеф Земли, область солнечного света для данной даты и времени в UTC (date):
[Код]
```
#скин физической карты Земли
ax.stock_img()
#добавление земной тени с параметрами date и прозрачностью в 0.2
ax.add_feature(cfeature.nightshade.Nightshade(date, alpha=0.2))
#отображение координатной сетки
ax.gridlines()
```
Приступим к тестам модуля. Как я упоминал ранее, рассмотрим 4 ситуации, взятые из личной практики, на которых разбираются теоретические вопросы. Будем придерживаться каноничной нумерации элементов Python.
Cartopy, но проверяй #0
Для начала убедимся, что точки, привязанные по координатам к Северному полушарию, не дублируются/не наслаиваются в Южном, т.е. что модуль адекватно обрабатывает и отображает данные по полушариям. Дабы точки были отличимы и более информативны, попробуем присвоить им сравнительную величину.
Отобразим данные через scatter (диаграмму рассеивания) – она хорошо применима в задачах представления объектов (точек) в значениях двух измерениях (абсциссе, ординате соответствуют долгота, широта объекта). В качестве данных, к примеру, возьмём что-нибудь связанное с деньгами — датасет старбаксов, а из него отобразим те магазины, широта которых больше нуля. В параметре «c» принято задавать массив значений, интерпретирующий какой-либо исследуемый количественный параметр (например, процент людей в каждом регионе, проголосовавших за кандидата в президенты, уровни землетрясений в различных эпицентрах, общая стоимость и объёмы продуктов) с помощью цветовой шкалы. Взглянем на датасет с помощью pandas. Число магазинов узнаем через shape – таких значений 25 600.
[Код]
```
starbucks_df = pd.read_csv(‘directory.csv’, delimiter=’,’)
print(starbucks_df)
print(starbucks_df.shape)
```К сожалению, сравнительного количественного показателя в датасете не нашлось. Ну и ладно. Сгенерируем посещаемость магазинов (человек за сутки) в нормированном виде – в отрезке от 0 до 1. Для генерации значений используем random.uniform с округлением до 2-х знаков после запятой. Количество таких значений равно 25 600 (числу магазинов).
Pop_vals = [round(x, 2) for x in np.random.uniform(0, 1, 25600)]Значения широт расположены на отрезке [-90;90], долгот: [-180;180]. Нас интересует всё по широте = (0;90].
[Код]
```
stbs_above_eq = starbucks_locations[starbucks_locations[«Latitude»] > 0] #в датафрейме элементы с широтой > 0Палитра цветов присваивается в параметре «cmap». Добавим списки в scatter:
[Код]
```
plt.scatter(lons, lats, c=pop_vals, cmap=’jet’, marker=’.’, lw=0.5, s=20, transform=ccrs.PlateCarree())
```Покажем точки на ортографической проекции под разным углом, дабы убедиться в корректности:
Как видно по итогу, Cartopy убрал полушарие, в котором не было точек. «Подрезал» и область около полярного полюса, т.к. Старбакс там пока не открыт. Также удалось сделать точки отличными между собой – по ним видно, что магазины в основном расположены в США, Европе и Восточной Азии. Поскольку данные о посещаемости рандомизированы, визуально определить области с наибольшей популярностью магазинов в данном случае затруднительно.
Cartopy, но проверяй #1
В задачах возникают случаи, когда нужно перейти от одной проекции к другой. Корректно ли перестроятся системы координат при смене проекции?
Для этого зададим координаты Иркутска, отобразим их на Plate-Carree, после сменим проекцию на Orthographic.
[Код]
```
#предварительно введены координаты в irk_lon, irk_lat
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=100.0, central_latitude=30.0))
plt.scatter(irk_lon, irk_lat, color='red', marker='.', s=100, transform=ccrs.Mercator())
plt.text(irk_lon+2, irk_lat+2, 'Irkutsk', color='blue', transform=ccrs.PlateCarree())
```
Стоит отметить transform — крайне важный параметр. Он задаёт, в какой системе координат первоначально определены данные, рекомендуется постоянно его указывать, т.к. crs может автоматически предположить «неожиданную» систему — результаты визуализации будут некорректными.
Для сравнения в качестве инструмента-эталона работы с координатами возьмём сервис Google Earth.
Для примера те же геоданные введём в сервисе Google Earth в качестве инструмента-эталона при работе с координатами на Земле. При наглядном сравнении точки расположены правдоподобно — «середина Земли», чуть западнее Байкала.
Cartopy, но проверяй #2
Возможно ли с помощью web-запросов и модулей перейти в более «локальные» масштабы и показать достопримечательность города?
Да, это оказалось возможным. В Cartopy входит модуль io (input/output), который позволяет получать изображения и карты: загружать, сохранять и извлекать карты в различных форматах данных. Классы io.img_tiles представляют собой интерфейсы к соответствующим веб-ресурсам (например, Google Maps и Stamen) для автоматической загрузки «тайлов» (фрагментов карты) — они нам и понадобятся. Используем класс io.img_tiles.OSM для работы с OpenStreetMap (популярной и свободной для использования географической картой мира). Помимо него необходимо подключить библиотеку requests для обращения по API, PIL — для получения изображения в tile‐формате. Реализация фрагментарно выглядит следующим образом.
[Код]
```
#функция для переформатирования веб-запрос от OSM для Cartopy
def get_prep_img(self, tile):
#url – заданный URL-адрес API карты улиц
url = self. _image_url(tile)
r = urllib.request.Request(url)
fh = urlopen(r)
#получаем изображение
img_data = io.BytesIO(fh.read())
#закрываем url
fh.close()
#открываем изображение с помощью PIL
img = PIL.Image.open(img_data)
#устанавливаем формат изображения
img = img.convert(self.desired_tile_form)
return img, self.tileextent(tile), 'lower'
#применяем Cartopy
cartopy.io.img_tiles.OSM.get_image = get_prep_img
osm_img = cartopy.io.img_tiles.OSM()
fig = plt.figure(figsize=(10,10))
#нужна проекция, которая использует систему отсчёта координат (CRS) StreetMap
ax = plt.axes(projection=osm_img.crs)
moscowgate_point = [lat, lon]
#эвристически подобранное отдаление от координат Московских ворот
zoom = 0.00075
extent = [center_pt[1]-zoom*2.0),center_pt[1]+(zoom*2.0),center_pt[0]-zoom,center_pt[0]+zoom]
# приближенный фрагмент карты получаем с помощью set_extent
ax.set_extent(extent)
# scale – параметр масштаба, подобранный эвристически
ax.add_image(osm_image, scale)
```Итоговое изображение будет следующим.
Если действовать по тому же принципу, но класс OSM (cartopy.io.img_tiles.OSM) изменить на QuadtreeTiles (cartopy.io.img_tiles.QuadtreeTiles), картинка станет более «реалистичной» — получится спутниковое изображение. Добавим ещё пару штрихов, которые к Cartopy не имеют отношения: поиграем с величиной масштаба, достопримечательности отметим крестиком, изобразим сетку и радиус-круг с геодезическим расстоянием в 1 км от центральной точки. Получившийся результат представлен ниже.
Cartopy, но проверяй #3
Помимо отображения объектов на карте, в различных задачах также порой нужно определить географические величины. Классический пример – найти расстояние между точками. Возник вопрос: способен ли Cartopy определять такое расстояние?
Задача: допустим, Николай Васильевич, матёрый аудитор и ярый любитель футбола, планирует посетить Чемпионат мира. При этом он тщательно следит за здоровьем и подсчитывает свои шаги ежедневно. Известно, что время = деньги и будет использован самолёт, но интересно и неизвестно, сколько шагов* «упустит» ревизор.
*Николай Васильевич соблюдает стандарты повсюду, поэтому имеет размер шага в эталонный метр.
С Cartopy нас ждёт промах – инструмента для обработки данных в пакете не обнаружил. Поиск проводил в документации актуальной версии модуля.
Решение:
Определим расстояние самостоятельно: первое, что приходит на ум – использовать для измерения евклидово расстояние. Пусть a и b – две точки, λ, и ϕ – соответствующие точке широта и долгота, a = (λ1, ϕ1), b = (λ2, ϕ2), тогда расстояние между ними:
Но для сферы/эллипсоида эвклидово расстояние недопустимо, оно приводит к существенным неточностям. В таком случае применяется формула Haversine, где кратчайший путь между двумя точками на сфере – это путь по дуге большего круга. 2 арксинуса необходимо умножить на радиус сферы r. В нашем случае это радиус Земли + высота над её поверхностью. Высота будет равна нулю, поскольку города расположены на поверхности. Итоговая формула примет следующий вид:
где r – радиус сферы (Земли).
Реализуем формулу в виде функции.
[Код]
```
def get_haversine_distance(point1, point2):
# формула тригонометрическая, поэтому десятичные градусы необходимо преобразовать в радианы
lon1, lat1 = [np.radians(x) for x in point1]
lon2, lat2 = [np.radians(x) for x in point2]
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
return 2 * EARTH_RADIUS * np.arcsin(np.sqrt(a))
irk_qat_dist = get_haversine_distance(irk_point, qat_point)
```Результат (irk_qat_dist) изобразим с помощью Cartopy — добавим в виде title.
Таким образом, Николай Васильевич упустит 5,314 млн. шагов.
Благодаря методу set_extent(lon1, lat1, lon2, lat2) можно показать только фрагмент карты, исходя из заданных координат. Для наглядности убрал отображение поверхности (stock_image).
Изобразили, проверили, теперь резюмируем.
Как использовать Cartopy — всё зависит от Вашей идеи и сферы поставленной задачи. Пример – финансовое учреждение. Если хотим показать процентное соотношение просроченных кредитов ко всем по территориальным областям, для этого хорошо подойдут pie’sы, если планируем изобразить новые структурные подразделения на карте и их функциональность в градации от красного к зелёному — смело можно брать scatter, потоки средств — streamplot, а contourf продемонстрирует области по уровню использования продуктов, а можно ещё реализовать bubble map, картограмму, карту связей и т.п.. Из моего опыта проекции Plate-Carree, Orthographic и диаграмма рассеяния scatter покрывают большинство задач по визуализации объектов на поверхности Земли, околоземном пространстве.
«Где карта, Билли?»
Из минусов стоит отметить возможные трудности при установке Cartopy — проблемы с совместимостью зависимостей, ошибки с wheel’ами и PEP. Могу посоветовать использовать conda и Linux/WSL. Также Cartopy пока предоставляет возможность только визуализировать данные, нет в нём обрабатывающих функций «из коробки», чтобы найти расстояние между объектами, площадь сектора и т.п. И, наконец, отсутствует интерактивный режим с картой — результат невозможно ни приблизить/отдалить, ни повращать, объекты на карте «некликабельны» — а хотелось бы, кликнув по ним, узнать дополнительную информацию, если визуализация размещается на веб-сайте.
Тем не менее, Cartopy — мощный инструмент, достойный внимания и применения, результаты которого пригодны для дальнейшего анализа и публикации, например, в научных работах. Пакет позволяет получать доступ к географическим структурам данных через такие форматы, как shapefile’ы и GeoJSON. В следующей публикации можем сравнить Cartopy с другими инструментами визуализации данных на картах (Basemap, Folium, у них свои карты и фишки), а можем рассмотреть методы компьютерного зрения и машинного обучения, применимые в одной реальной исследовательской задаче – входными данными послужат изображения, полученные как раз с помощью Cartopy. Ваши пожелания по выбору темы оставляйте в комментариях. Там же готов обсудить вопросы, рекомендации, касаемые данного материала.
D.S. Для интерактива прикрепляю гифку: здесь с интервалом в 5 минут измеряется X. Смело оставляйте предположения, что скрыто под X. Первый угадавший получит приз — нет, не карту с деньгами, а полный код этой публикации.
